Universal structures for adaptation in biochemical reaction networks

At the molecular level, the evolution of life is driven by the generation and diversification of adaptation mechanisms. A universal description of adaptation-capable chemical reaction network (CRN) structures has remained elusive until now, since currently-known criteria for adaptation apply only to a tiny subset of possible CRNs. Here we identify the definitive structural requirements that characterize all adaptation-capable collections of interacting molecules, however large or complex. We show that these network structures implement a form of integral control in which multiple independent integrals can collaborate to confer the capacity for adaptation on specific molecules. Using an algebraic algorithm informed by these findings, we demonstrate the existence of embedded integrals in a variety of biologically important CRNs that have eluded previous methods, and for which adaptation has been observed experimentally. This definitive picture of biological adaptation at the level of intermolecular interactions represents a blueprint for adaptation-capable signaling networks across all domains of life, and for the design of synthetic biosystems.

where g(x, y) = 0 and y is a non-RPA capable variable which forms the kinetic pair to x. It is unclear why g(x, y) can only be a function of one additional variable, apart from the output variable. The proof given on page 12 in the Supplement does not satisfactorily explain why the ideal I f ∩ R[x] will contain polynomials in x j and x m that are not in I f ∩ I p .
As an example consider the following network: Furthermore assume that species X 2 , . . . , X 5 participate is several reactions that do not involve X 1 but can be catalysed by it. In this scenario the RPA polynomial would be ρ = x 2 x 3 ( c 1 c 2 − x 1 ). Hence the function g depends on two non-RPA variables x 2 and x 3 . Please explain how this example is consistent with the form of g(x, y) stated above.
Moreover, it is not immediately clear how existence of multiple opposer/balancer modules translates into existence of existence of corresponding RPA polynomials. This has not been explained in sufficient detail.
• Connection to existing works: In [4] the authors consider RPA systems that are maximally robust and find simple linear-algebraic structural conditions that characterise this property in both deterministic and stochastic settings. How do the results in this paper connect to the results in [4]?
Importantly, our new and unexpected findings on the relationship between the deficiency of a CRN and its overarching topological structure allow us to make a new and crucial connection to our previous work on RPA at the network macroscale [1], thereby enabling us to identify, definitively, the one-to-one correspondence between the subsidiary linear control problems and key topological features of the CRN. The fact that we can connect these new results to our previous general (macroscale) topological solution to the RPA problem is, without question, extremely powerful. But it is entirely untrue that we simply impose our previous results on RPA topologies at the network macroscale onto the microscale of intermolecular interactions in a CRN, leading to a result of "incremental novelty" as the reviewer goes on to suggest. In fact, no previously published work to our knowledge has been able to identify the general relationship between deficiency and linear complex invariants in RPA-capable networks. Cappelletti et al. [2], for instance, do note that an "RPA polynomial" (as we call it here) can exist in the rowspan of a CRN for deficiencies greater than one -if the RPA polynomial in question happens to be in the rowspan of the CRN, which is not generally the case -BUT fail to identify how any added deficiency (beyond unity) allows this crucially important "RPA polynomial" to reside in the rowspan. As a consequence, even the very simplest RPA-capable CRNs whose RPA polynomial does not reside in the rowspan (e.g. the toy model discussed by those authors, that we analyse as Example 1 in our SI (Section S3.1)) simply elude the Cappelletti et al. [2] approach. No previous work as been able to characterise, in complete generality, the nonlinear coordinate changes that are needed to reveal the capacity for RPA in CRNs in an entirely universal manner.
We are truly grateful for the reviewer's fresh set of eyes on our work, and have carefully revised our manuscript (and SI) to ensure that the novelty and scientific importance of our findings have been clearly and accurately communicated. We acknowledge, in particular, that the preamble to our Supplementary Information (i.e. the Supplementary Text prior to Section S1) could have been organised much better, with the aims and novel contributions of our work expressed much more explicitly. We acknowledge that our original submission mentioned the development of a definitive 'RPA test' first in this section, before describing the overarching goal of demonstrating the truly universal principles of RPA in CRNs, and the relationship of this universal description of RPA to integral control. This section has now been carefully revised to read as follows: The overarching goal of this work is to provide a comprehensive and universal description of all possible Chemical Reaction Network (CRN) structures that can implement the keystone biological function known as Robust Perfect Adaptation (RPA), and in so doing, to propose a definitive algebraic condition that captures presence of the RPA property in any CRN that possesses it. We also recognise here that all forms of RPA, including the special case known as Absolute Concentration Robustness (ACR), must implement some form of integral control. Our universal characterisation of RPA at the level of intermolecular interactions thereby provides an algorithmic method not only for detecting the RPA property, but also for identifying the hidden integral controllers, as well as the 'setpoints' of any RPAexhibiting molecules.
As we show in the pages to follow, providing such an all-encompassing description of RPA at the microscale-level of biological networks -that is, in CRNs, which account for the mathematical graph structure of intermolecular interactions -has required us to reconcile a wide range of disparate mathematical viewpoints, ranging from Chemical Reaction Network Theory (CRNT), control theory and the internal model principle, and algebraic geometry. The essential structure of our mathematical development is as follows: 1. In Section S1 we identify a universal 'kinetic pairing' principle, which captures the essence of RPA in all CRNs that possess the property, and which is codified by our Two-Variable Kinetic Pairing Theorem (Theorem 1, Section S1.4). This theorem guarantees that, for all RPA capable CRNs, the geometric projection of the CRN's steady-state locus onto a state space defined by just two variables is captured by a distinguished algebraic invariant that we call an RPA polynomial. In this context the term 'variable' has a particular meaning, that we define carefully in Section S1. 4.
In this Section we also introduce readers to the fundamentally important notion of algebraically independent subnetworks -a simple yet powerful concept introduced by Martin Feinberg (see Appendix A.6 in [9] for a detailed overview), whose implications for the decomposition of CRNs into independent subnetworks, and ultimately topological modules, that are compatible with RPA are recognised for the first time in the present work.
This key concept provides a vital ingredient which contributes to a definitive characterisation of RPA-capable CRNs of arbitrarily high deficiency, far beyond the simple deficiency-one ACR-capable CRNs covered by the Shinar-Feinberg theorem [7].
2. In Section S2 we briefly discuss the topological implications of Theorem 1 and, for later reference, make a preliminary connection with previous work on the topological structures that are known to hold, at the network macroscale, for all RPA-capable networks [1]. Importantly, the universal topological solution to the RPA problem [1] was developed considering only a single external input or disturbance, and did not account for the possibility of no external inputs, with alterations in total abundances of the interacting elements (as specified by the initial conditions of the system) constituting the only possible perturbation to the system. More importantly, and in stark contrast to the present work, the universal macroscale solution to the RPA problem [1] gives no concrete information as to the CRN structures (at the network microscale) that could engender RPA or satisfy its strict topological requirements.
3. Since most RPA-capable CRNs identified in prior studies identify an RPA polynomial (as we call it here) from only linear coordinate changes (see, for example, Cappelletti et al. [2]), in Section S3 we consider the algebraic properties of CRNs that allow the RPA polynomial (as expressed in Theorem 1) to reside in the rowspan of the CRNs chemical reaction rates (i.e. given by an ℝ-linear combination of the rate equations). For this purpose, we analyse two carefully chosen CRN examples in detail, and demonstrate the existence of several subsidiary polynomial invariants in the rowspan for each case, and note the relationship of these polynomial invariants to the overarching topology of the system. We demonstrate via deficiency-preserving transformations of the respective CRN structures that these subsidiary polynomial invariants combine to yield the all-important RPA polynomial within the rowspan exactly when the invariants are stoichiometrically dependent. If the invariants are stoichiometrically independent, on the other hand, a concatenating monomial is required to reconcile the invariants, thereby introducing a necessary nonlinearity to the transformation that produces the RPA polynomial. The detailed analysis of these simple examples provides the reader with an accessible overview of the general principles that characterise RPA-capable CRN structures before we develop these principles more rigorously and in complete generality in the remaining Sections. 4. In Section S4 we demonstrate, through an analysis of deficiency-increasing reactions within independent CRN subsets, that the coordinate change required to extract the all-important RPA polynomial from the CRN rate equations is always almost linear, in the sense that the RPA polynomial can always be decomposed into a collection of 'complex linear invariants'. From a control theory viewpoint, these observations lead to the previously unrecognised conclusion that the 'internal model' for any RPA-capable CRN can always be decomposed into a collection of linear integral controllers, in addition to reconciling the processing of biochemical information at the network microscale (i.e. via the intricate intermolecular interactions that comprise a CRN) with the overarching topological features at the network macroscale [1]. 5. In Section S5, we provide the reader with extensively annotated code in the open-source software Singular (www.singular.uni-kl.de) to test the RPA capacity of a range of illustrative examples. This code can be adapted readily to the analysis of any CRN. If RPA does obtain for a particular CRN under investigation, an important consequence of the special almost-linear structure of the transformation required to identify the RPA polynomial is that the computational demands of this algorithm do not differ substantially from those of Gaussian elimination (a polynomial-time algorithm).
We clarify the nature of our scientific advance in further detail as we address each of the reviewer's individual points in turn.
• Lack of sufficient novelty: -The general principles for RPA have been developed by (slightly) extending previous results by the authors for deficiency zero and deficiency one CRNs. These low-deficiency networks are related to `balancer' and `opposer' modules which are known to form a topological basis for all RPA networks [1]. This connection is explained with two illustrative examples and general principles are formulated by relying on the theory in [1]. Hence, given the earlier results in [1], the incremental novelty in this paper is not significant.
Our response: We wish to respectfully emphasize, first and foremost, that these completely general principles by which RPA is always realised at the level of intermolecular interactions in a CRN have by no means "been developed by (slightly) extending previous results (of ours) for deficiency zero and deficiency one CRNs". Our framework is certainly not limited to "low-deficiency networks", as the reviewer claims. The general principles, which we succinctly summarise for a general audience via two particularly simple illustrative examples, in order to provide as accessible an overview of the main ideas as possible at the outset, have most assuredly not been "formulated by relying on the theory in [1]".
Our previous work [1], considers the broad macroscale topological structures that are compatible with an RPA response, and delineates the full solution space to the RPA problem from this viewpoint. This approach considered an abstract network, with a single 'input stimulus', or disturbance, delivered to a distinguished input node. The network was considered to exhibit RPA exactly when a distinguished output node, not necessarily distinct from the input node, could always return to a single, fixed value at steady-state (the 'setpoint'), independently of the magnitude of the input stimulus, and independently of the choice of parameter values. Using a topological method, where the flow of biochemical information could be partitioned into the subsets of a topology that accounted for the algebraic requirements of RPA, an exhaustive description of all possible RPA-compatible network designs, induced relative to the chosen input-output node pair, could be identified. This approach definitively established that all RPA-capable networks of this type must be decomposable into two types of well-defined topological structures, or 'modules' -Opposer modules and Balancer modules -and concluded that certain key 'computational nodes', corresponding to well-defined constraints on reaction kinetics, must be embedded into these modules. In particular, Balancer modules require one or more balancer nodes, exhibiting 'balancer kinetics' along with a single connector node, exhibiting 'connector kinetics'; Opposer modules, by contrast, require one or more opposer nodes, exhibiting 'opposer kinetics', organised (topologically) into structures known as 'opposing sets'.
But these conclusions, by themselves, provide absolutely no information on how these 'reaction constraints' might be realised, in general, by the intricate interactions among collections of molecules. Nor has it been clear until now how any of the known RPA-capable (or ACR-capable) collections of chemical reactions (CRNs) encode any of the topological requirements for RPA that we had identified in [1]. Indeed, in our earlier paper [1], our illustrative examples employed simple 'invented' functional forms (disconnected from CRN structures) which encapsulated the principles of opposer kinetics, or balancer/connector kinetics, in the most straightforward way, since the focus of that earlier work was the overarching network design principles that must hold at the network macroscale. An opposer node , regulated by some node , for instance, could be assigned a kinetics of the form / = − , which is zero-order in . But this contrived reaction form could certainly not be induced by a CRN for a variety of mathematical reasons, including the fact that (by the so-called 'Hungarian Lemma'), all terms preceded by a negative sign must include the species that is the subject of the reaction equationin this case, . Indeed it has been entirely unclear, until now, how the necessary conditions for RPA could be realised in general by 'real' chemical reactions involving molecules that interact in specific and potentially intricate ways according to a mathematical graph structure (i.e. a CRN). Previous studies (e.g. Shinar and Feinberg [7], Cappelletti et al. [2], Karp et al. [8], etc.) have only provided partial answers for particular special cases.
Of course, the antithetic integral control motif is a simple and well-known reaction structure that can serve as an 'opposer node', and can encode opposer kinetics (and thereby implement integral control) in a fairly obvious and straightforward manner via a particularly simple linear transformation. But this mechanism by no means explains the RPA (ACR) property observed in the Shinar-Feinberg deficiency-two model of the EnvZ-OmpR network, for example. In fact, until now, it has not been at all clear how the EnvZ-OmpR motif, or virtually any other known example of an ACR/RPA-capable CRN, might relate to the topological principles we identify in [1]. That the completely general description of RPA-permissive CRN graph structures, as we introduce in the present paper, should be able to map in such a clear, direct and well-defined way with our earlier universal description of RPA at the network macroscale is both remarkable and immensely satisfying. In particular, the revelation that balancer 'nodes', opposer 'nodes' and connector 'nodes' must necessarily be linear transformations of the CRN's molecular concentrations and that, combining these invariants nonlinearly in the very specific manner we describe here, gives rise to a completely general description of all possible RPA-capable CRNs, is entirely unexpected. A priori, there was no reason to think that there would necessarily be such a straightforward way to concretely reconcile the CRN microscale with established design principles of the network macroscale [1] on the basis of any previous work. Only in retrospect does this seem inevitable! In any case, no previous study has made the well-defined connections we make here between CRN deficiency, network topology (encoding the flow of biochemical information through the graph structure of a CRN) and the internal model principle, in addition to the general results pertaining to the decomposition of the requisite nonlinear map (that extracts a suitable internal model) into a collection of linear maps, governed by the fundamental algebro-geometric requirements of RPA (and hence, ACR). As a consequence, no previous study has been able to provide a definitive description of all possible RPA-capable (and, more specifically, ACRcapable) CRNs, nor a well-defined algebraic condition (and associated algorithmic test) that captures the presence of the RPA property in any CRN that possesses it. Although there have certainly been previous ad-hoc attempts (e.g. in [5]) to demonstrate ACR/RPA capacity in a specific CRN via computation of a suitable Gröbner basis, these ad-hoc attempts, of themselves, are not generalisable (as we explain in greater detail later, in response to the reviewer's more specific queries regarding prior work). Moreover, there was no reason to think that this approach would be computationally feasible in general prior to our demonstration of the universal decomposition of the underlying nonlinear coordinate change into a connected sequence of linear maps. (Again, this is what we mean when we claim that the nonlinear map is 'almost linear'). Of even greater computational concern is that fact that, prior to our proof that the RPA polynomial (as we call it here) is always obtained via geometric projection onto exactly two variables (noting that 'variable' does not necessarily mean 'species'!), prior ad-hoc attempts to compute RPArelevant Gröbner bases have been forced to use the most computationally expensive monomial ordering (i.e. the lexicographic ordering), since it was not known at the outset which, or how many, variables had to be eliminated. These former potential computational challenges have all been robustly resolved by the present work.
Of course, most of the RPA-capable CRNs that have been identified at the present time exhibit a special form of RPA known as Absolute Concentration Robustness (ACR). Unlike the abstract networks we consider in [1], there is no external 'input stimulus' for these networks, whose steady-states can only be altered by varying the total abundances of the various interacting molecules, as specified by the initial conditions for the reactions. For example, the Shinar-Feinberg models of the EnvZ-OmpR motif, comprising both deficiency-one and deficiency-two versions, are known to exhibit ACR in the phosphorylated OmpR moiety, yet no previous study of this motif (including the recent analysis of the deficiency-one version by Cappelletti et al. [2]), has succeeding in connecting the ACR capacity of this particular molecule to its topological balancer structure (and certainly not in a way that allows a general connection to be made to the set of all possible RPA/ACR-capable CRNs). As we show in great detail in our newly revised manuscript, a balancer structure arises from the fact that, for this CRN, the one unit of deficiency encodes two parallel pathways that regulate the phosphorylation of OmpR, whereas the two units of deficiency (in the deficiency-two version) encodes three such feedforward pathways that regulate the same phosphorylation cycle. Indeed, it is the deep connection we identify between deficiency and the flow of biochemical information (topology), along with the relationship to modularity of the decomposition of RPA-capable CRNs into algebraically-independent subnetworks, that feature among the major conceptual advances of the present work. In our original submission, we had devoted the entirety of sections S4.2.1 and S4.2.2 in our SI to a succinct consideration of these ideas. In our revised submission, we have extensively revised the preliminary sections of our SI, adding in significant additional discussion on the significance of deficiency (a key CRNT concept) in Section S1.2, and have added in an entire new section (Section S1.3), which explains the concept of algebraically independent subnetworks of a CRN in great detail. (This is in contrast to our original submission, in which our far briefer discussion of algebraically independent subnetworks was delayed until Section S4.2.1). In our revised SI, the crucial connections between deficiency and topology are thereby made clear and explicit from the outset. We also briefly analyse the example of the Cappelletti et al. [2] toy model from the point of view of a decomposition into algebraically-independent subnetworks, along with a careful consideration of the partitioning of deficiency and rank into these subnetworks, so that the topological discussion of this (and other) networks in later Sections, as well as the general relationship of CRN topologies to deficiency and independent subnetworks, can be more readily appreciated by the general reader.
Importantly, ACR-capable CRNs, despite not being driven by exogenous inputs, nor perturbed by exogenous disturbances delivered to specific molecules, are still subject to the topological requirements that unify the broader class of RPA-capable CRNs which generally do involve external inputs -a deeply important result which fundamentally governs how such CRNs respond to perturbations by other (new) classes of disturbances (e.g. enzyme inhibitors and other small molecules), or new exogenous inputs/disturbances distinct from those represented in the original CRN. Since the detailed reaction structures governing complex cellular signalling networks (e.g. signal transduction networks in Metazoa) remain mostly unknown, having access to such a complete and universal description of RPA-capable CRN structures in the abstract is immensely powerful, particularly for considering new possibilities for molecular-targeted therapies, for instance, and for considering evolutionary alterations to CRNs.
We do acknowledge that our earlier work [1] foreshadowed the principle of a topologically distributed internal model, where one internal model could 'feed into' another topologically related one. Nevertheless, no prior work has been able to reconcile the mathematical properties of dynamical systems induced by CRNs with these previously-identified principles, and our current work highlights that making the fundamental connection between macroscale (e.g. from [1]) and CRN microscale poses significant mathematical challenges. Our new study clearly delineates the full suite of new mathematical ideas required -on the relationship of subsidiary linear invariants (polynomials) to the algebraic elimination problem; on the relationship of stoichiometric independence (which is itself a function of the algebraic structure of the vertices of the CRN graph) to the inability of polynomial invariants to 'pass' within the rowspan; and on the requirement for 'concatenating monomials' and the relationship of these to the topological structure of the CRN (for which we identify a previously unrecognised connection with CRN deficiency, as we explain in far greater detail in our revised submission).
We truly value the reviewer's feedback, which has helped us appreciate where our original submission failed to compellingly communicate the novel aspects of this work. We are deeply grateful for the reviewer's time and effort in re-considering our extensively revised work, and are thankful for his/her support.
-Many researchers have exploited Gröbner basis and ideal construction for obtaining steady-states in terms of parameters and in particular for studying ACR. For example, see [6,5,3] and references therein. It seems that the authors are unaware of these works and in order the assess the novelty of the methods in this paper, these existing works must be cited and discussed.
Our response: We thank the reviewer for highlighting these references, which use Gröbner bases to elucidate various aspects of CRNs and even consider ACR in one instance [5]. We value the reviewer's suggestion that the novelty of our new work could be made clearer by citing and discussing earlier work employing Gröbner bases and ideal construction in the context of chemical reaction networks.
We certainly acknowledge that the elucidation of polynomial ideals via Gröbner basis computations is well-known to have wide-ranging applications to problems in science and engineering [10]. We do not claim that the use of Gröbner bases, in and of itself, to elucidate some property of a polynomial system is novel.
There are actually exceedingly few published papers in the literature currently that employ Gröbner basis methods, or 'ideal construction' to study RPA/ACR specifically, other than the article by Perez Millan et al. [5]. Those authors compute a Gröbner basis in an ad-hoc manner for the deficiency-two EnvZ-OmpR motif, and use this result to demonstrate the capacity for ACR in one of the species. But the chief focus of that paper [5], and of the Craciun et al. paper [3], is not ACR, or RPA, but CRNs that exhibit toric steady-states. As our study makes clear, RPA-capable CRNs do not, in general, give toric ideals, or have toric steady-states. Although the EnvZ-OmpR motif analysed by Perez-Millan et al. [5], along with several of the illustrative examples we consider in our paper are characterised by toric ideals, it is clear that this could not possibly be a general property of RPA-capable CRNs for a number of reasons, including the fact that the pairing function ( , ) (see our Theorem 1) can be a multi-term polynomial, and not simply a power product of species (monomial). In addition, even in cases where the RPA polynomial itself is a binomial, this condition does not guarantee that the steady-state ideal has a (standard) generating set composed entirely of binomials. A simple example of an RPA-capable CRN that is not generated by binomials, and is therefore not characterised by a toric ideal, is given in Figure 2 (previously Figure 1) of our main paper. (See, in particular, SI S5.2, where we compute a number of relevant Gröbner bases for this particular CRN, thereby demonstrating that the ideal corresponding to its chemical reaction rates cannot be generated by binomials). Thus, while the class of RPA-capable CRNs has a non-empty intersection with the class of toric dynamical systems, the study of toric ideals, unfortunately, cannot capture the fundamental essence of the RPA property in CRNs.
In our extensively revised SI, we provide a detailed commentary on the computational demands of Gröbner basis-computing algorithms from a complexity standpoint (see Remark 3 following Definition 3 in Section S1.5), and on the simplifying structure afforded by RPA-capable CRNs (owing to the 'almost-linear' coordinate change that we demonstrate can always extract the 'RPA polynomial' from the CRN rate equations). At the conclusion of those explanations, we explain that " … the ACR-capacity of the deficiency-two EnvZ-OmpR motif was demonstrated via computation of a Gröbner basis, in an ad-hoc manner, by Perez Millan et al. [5]. But without knowing which, or how many, variables generally characterise the RPA-encoding elimination ideal, and without establishing the special 'almost linear' structure of the mathematical transformation required to identify the RPA polynomial in any RPA-capable CRN, this approach does not generalise, and cannot be applied to the systematic analysis of large and complicated CRNs (e.g. in the context of metabolism [11] for which candidate RPA-capable molecules have not already been posited from experimental evidence)." Note that this is not a criticism of the Perez-Millan et al. [5] paper in any way. Again, the focus of the Perez Millan et al paper is not ACR-capable networks, so there is nothing amiss with their ad-hoc approach to the particular CRN they consider, which happens to be characterised by a toric ideal (which is the focus of their paper); they do not claim to present a general method that will necessarily be valid (or computationally feasible) for any ACR-capable CRN. They simply show that the particular CRN they consider does exhibit ACR, and the particular Gröbner basis they compute confirms this. The Sadeghimanesh and Feliu paper [6], on the other hand, considers the computational demands of Gröbner basis algorithms for the analysis of CRNs, and proposes a transformation of CRNs containing 'intermediate species' to yield a core network from which these intermediates are absent. This approach allows a more computationally efficient monomial order to be employed in Gröbner basis computations, where only the intermediates are subject to the computationallyexpensive lexicographic ordering, and with faster orderings (e.g. grevlex -graded reverse lexicographic) applied to non-intermediate species. Thus, what these authors effectively achieve is the development of a particular 'block ordering' on monomials, and which can potentially speed up Gröbner basis computations quite significantly in comparison with a 'complete' lexicographic monomial ordering. This approach can be employed, in principle, for any CRN containing intermediate species (as defined in their paper [6]) -not just RPA-capable CRNs.
Our paper, by contrast, focusses squarely on RPA-capable CRNs, and aims to discover their general properties. For these CRNs, in contrast to the much broader class of CRNs considered by Sadeghimanesh and Feliu [6], the computational demands confronting any Gröbner basis computations are comparatively very light. There are two reasons for this: (1) If the CRN is indeed RPA-capable, then Theorem 1 guarantees that the RPA polynomial can be obtained via elimination of all but two variables. As a consequence (and as we point out explicitly in Section S5 where we instruct the reader on the implementation of these algorithms in Singular, using a variety of illustrative examples), we can always use a 'block ordering' on monomials, where the two chosen projection variables are ordered 'lower' than the remaining variables. Having established this, the computationally expensive lexicographic ordering is not required: both of the 'blocks' in this special elimination order can now accommodate a fast monomial ordering (e.g. grevlex, sometimes referred to as degrevlex, which uses the syntax 'dp' in Singular). But there is a second, and more important reason: (2) The transformation that reveals the RPA polynomial is 'almost linear' in the sense that only ℝ-linear combinations of the rows of the system are required to identified the 'fundamental' (subsidiary) polynomial invariants (which are subsequently combined nonlinearly via concatenating monomials). For this reason, the syzygy polynomials (S-polynomials) that are calculated as part of Buchberger's algorithm are identical to the linear row combinations that are computed during Gaussian elimination. But note that this computational 'streamlining' is linked inextricably to the special CRN structures that implement RPA, and establishing this key property is one of the highly novel contributions of our work. For more general CRNs, computing Gröbner bases may be computationally infeasible (especially if imposing a global lexicographic ordering at the outset, as one is often forced to do without any prior analytical guideposts); from this point of view, the approach offered by Sadeghimanesh and Feliu [6] is both interesting and potentially helpful.
We also emphasize that the existence of an algorithm that can test for the RPAcapacity in a particular CRN does not, in and of itself, reveal the definitive properties of all possible RPA-capable CRNs. Our paper is not "about" Gröbner bases as such. That we can, indeed, employ Gröbner basis-computing algorithms (e.g. Buchberger's algorithm) to test for RPA capacity, in general, is a consequence of the special properties of RPA-capable CRNs that we identify for the first time in this new paper.
-More importantly, the central analysis in this paper crucially depends on the result for deficiency one networks which is stated as Theorem 3 in the Supplement. However, this result is not new, as it follows from Theorem D.1 in the Supplement of [2] (see the proof in [2]).

Our response:
We have carefully scrutinized the proof to Theorem D.1 in the Supplement of [2], and acknowledge that the reviewer is entirely correct to note that our Theorem 3 follows directly from Theorem D.1 in [2], and that our essential mathematical argument is virtually identical to that in [2]. We have now updated our preamble to Theorem 3 (Section S4.2) to emphasize to the reader that this alternative version of Shinar-Feinberg's theorem draws on arguments developed by Cappelletti et al. [2], and provide a reference to Theorem D.1 in [2]. In the interests of a more comprehensive and scholarly exposition of the key ideas (for later use in the remaining sections of the SI), we have also taken this opportunity to add in a second version of the theorem (Theorem 4), applicable to deficiency-zero CRNs, which shows that binomials involving any pair of terminal complexes from the same terminal SCC must also reside in the rowspan of the CRN reaction equations. Following these two theorems, we also add a remark, again citing Cappelletti et al. [2] Proposition C.5, that the same mathematical arguments (as employed in Theorems 3 and 4) can readily be extended to demonstrate that binomials involving either a terminal complex and a non-terminal complex, or two terminal complexes from different SCCs, can never reside in the rowspan of a CRN.
We sincerely thank the reviewer for his/her careful attention to this key technical detail, and for ensuring that all mathematical concepts presented in our work are correctly attributed. We feel that the presentation of the mathematical ideas in this section are now much more comprehensive and scholarly, and the additional references to the Cappelletti et al. [2] paper are immensely valuable. We would like to respectfully emphasize, however, that the novelty of our present paper lies in not so much in our Theorem 3 (or D.1 in [2]) but in demonstrating how this result on deficiency-one (and deficiency-zero) networks may be extended to RPA-capable networks of any deficiency. The remainder of the SI, following Theorems 3 and 4, is devoted to precisely this goal. (The generalisations the arbitrary deficiency are also much more clearly and thoroughly communicated in our revised manuscript, with the help of the comprehensive explanations in the earlier sections of the SI (particularly Sections S1.2 and S1.3) on how deficiency governs the flow of biochemical information within a CRN (and ultimately relates to topology), as well as the decomposition of a CRN into algebraically independent subnetworks). We now include the following text at the beginning of Section S1 in our Supplementary Information:

Definition of RPA
In this study we consider Robust Perfect Adaptation (RPA) from the most general possible viewpoint. In particular, a CRN exhibits RPA in the concentration, , of some molecule exactly when maintains a constant steady-state value, (the molecule's 'setpoint'), for all steady-states of the system. The setpoint, , is a function of some collection of CRN parameters. Moreover, the CRN exhibits RPA (in ) in response to any perturbation or disturbance that does not feature in the functional form of the setpoint, .
With this very broad definition of RPA in mind, we recognise that there are many possible types of perturbations/disturbances that could alter the steady-state of the system: 1. One or more 'external' inputs. A disturbance of this type could arise in the form of some input molecule, , whose concentration is given by a step function generated by an 'exosystem' ( / = 0, not included among the CRN reaction rates). 2. Alterations in total expression levels (abundances) of the interacting molecules. If this is the only type of perturbation possible, then the CRN is considered a 'conservative network', corresponding to a 'closed' system (see [9], Chapter 4). In this scenario, the CRN captures the interconversion of molecules among a variety of possible forms. The dynamical system is thereby constrained to evolve within a particular stoichiometric compatibility class (an affine space parallel to the CRN's stoichiometric subspace -see Section S1.2) determined by the initial conditions. The orthogonal complement to the stoichiometric subspace determines the conservation relations for the CRN, which specify the constant expression levels of the various interacting molecules. Where a molecule can exhibit a fixed concentration across all steady-states when the CRN is subject to perturbations to these total expression levels, this type of RPA is referred to as Absolute Concentration Robustness (ACR, see Chapter 9 in [9]). 3. A disturbance that is encoded in a CRN parameter such as a production rate, , of some molecule, , as reflected in a CRN reaction of the form → . Note that CRN parameters determined by the intrinsic chemical properties of the interacting molecules (e.g. association/dissociation constants, catalytic constants, etc.) are very stable and not readily perturbed, except via mutation of the interacting molecules.
In the pursuit of greatest generality, we make no assumptions a priori as to which type of network perturbation might affect a CRN, nor do we impose any restrictions on which (or how many) parameters determine the setpoint. Of course, if the setpoint for a putative RPA-capable molecule involves a parameter that can be perturbed (see case 3. above) then it is clear that the molecule cannot exhibit RPA in response to that particular disturbance.
We can appreciate that the reviewer might have a different perspective on RPA from us, and hope that he/she is willing to accept our viewpoint and mathematical goals.
Regarding the reviewer's comment that "the set-point c (being) a rational function of the biochemical rate constants … would be the case for any steady-state" -we agree completely with the reviewer on this point. All variables of the CRN would certainly achieve some steady-state that is a rational function of parameters; but not all variables would achieve the same rational function of parameters across all possible steady-states of the system. By noting that the setpoint is a rational function of parameters, we do not intend to convey that the steady-states of other variables are not rational functions of parameters. We are simply pointing out the mathematical nature of the setpoint (which can ultimately be computed automatically via the algorithm we propose) for the benefit of the general reader, to whom such properties may not be obvious. We did not intend this to be a major point.
• Algorithm may not terminate: The Gröbner basis algorithm to find the RPA polynomial may not terminate. It is mentioned that failure to terminate for a chemical reaction network (CRN) is a prima-facie evidence that the CRN does not exhibit RPA. However this is not mathematically shown. In particular, no proof is given showing that non-termination implies no RPA. Moreover, non-termination of a method is not something that can be ascertained with any certainty by running the code. Accordingly, the authors cannot claim that they have a full 'characterisation' of the RPA property.
Our response: We thank the reviewer for this comment and for giving us the opportunity to provide greater clarity on these points in our revised manuscript.
First, we respectfully point out that Buchberger's algorithm, which computes Gröbner bases, always terminates. This is because all ideals of a polynomial ring (or any Noetherian ring, for that matter) are finitely generated -a result of central importance in algebraic geometry, formalized by the Hilbert Basis Theorem. As a consequence, we are assured that Buchberger's algorithm will terminate in a finite number of steps. The problem, of course, is that the finite number of steps required could, and often will, be unmanageably large! The computational demands of Buchberger's algorithm are well known, and we do acknowledge in our paper (in both our original and revised submissions) that the problem of determining whether a general collection of polynomials constitutes a Gröbner basis (also known as the Gröbner basis detection (GBD) problem) is NP-Hard. The computational demands of Buchberger's algorithm are clear from an examination of the algorithm itself; the demands are even greater when a lexicographic monomial ordering is used since this necessarily implements a computationally expensive polynomial division process in comparison with other monomial orders. The complexity class of the GBD has been well known since the 1990s (due mainly to the work of Bernd Sturmfels and others), and has been widely discussed in the literature.
We respectfully point out that our study does, in fact, make clear why RPA-capable CRNs have a special property which guarantees that Buchberger's algorithm will terminate in polynomial time if the CRN does indeed exhibit RPA. Indeed, we demonstrate that the mathematical transformation required to obtain the RPA polynomial from the reaction rate equations is almost linear in the sense that only ℝ-linear combinations of these equations are required to produce the 'subsidiary' polynomial invariants which are then 'concatenated' by monomials. And as long as only ℝ-linear combinations of the polynomials are required, the syzygy polynomials ("S-polynomials") that are computed at each step of Buchberger's algorithm are precisely the linear row combinations computed during Gaussian elimination (a polynomial-time algorithm). Indeed, for the simple CRNs considered by Cappelletti et al. [2], for instance, where the RPA polynomials of the CRNs were in the rowspan (no concatenating monomials required), then Buchberger's algorithm is, in fact, Gaussian elimination. Those authors present an efficient method to calculate which particular linear combination of the reaction rates is required to 'achieve' this RPA polynomial. But it is clear that for any CRN with an RPA polynomial in the rowspan of its reaction equations, the RPA polynomial itself can be obtained algorithmically by Gaussian elimination. In particular, if the equations are expressed in the matrix form dX/dt = Ax, where x is the vector of monomials, and the two monomials of the candidate RPA polynomial are listed last (at the bottom of the vector), then performing Gaussian elimination on A will result in an echelon form in which the bottom non-zero row has exactly two entries (corresponding to the two coefficients of the constituent terms of the RPA polynomial). Which linear combinations are required to achieve the rows of this echelon form can be stored during the execution of Gaussian elimination (producing an LU-decomposition, for example).
This process readily extends to the more general (nonlinear, but 'almost linear') case. The row combinations (corresponding to construction of S-polyomials) that are calculated during the execution of Buchberger's algorithm may be stored, and retrieved subsequently using the 'LIFT' command in Singular. We demonstrate this in the detailed annotations to our code in Section S5. Where only ℝ-linear combinations of the CRN equations are required to identify the RPA polynomial (as is for the cases considered by Cappelletti et al. [2]), the coefficients retrieved by LIFT will contain only 'constants' (polynomial functions of the CRN parameters). If the RPA polynomial is not in the rowspan, the row combinations retrieved by LIFT will contain (only) the concatenating monomials required to combine adjacent 'complex linear invariants' (subsidiary polynomials corresponding to topological features of the underlying flow of biochemical information in the CRN), in addition to the constants involving CRN parameters.
In reflecting on the reviewer's feedback, we have made extensive revisions to our paper to ensure that the technicalities surrounding the computational complexity of Gröbner basis computations for RPA-capable networks are completely clear and explicit. We now make extensive clarifications on these points throughout the SI. In Remark 3 following our proof to Theorem 1, for example, we explain that "…With a well-defined algebraic elimination problem now at hand, with two candidate variables selected for the geometric projection, the existence of an RPA polynomial as the sole generator of the relevant elimination ideal can always be tested in principle via the computation of a suitable Gröbner basis [10]. Unfortunately, the problem of computing Gröbner bases (e.g. via Buchberger's algorithm) for general polynomial systems is well known to be NP-Hard [refs], so without any special simplifying structure to the collection of polynomials, we can have no assurance that any Gröbner basis-computing algorithm will terminate on a practical time-frame. However, we will show in the Sections to follow that RPAcapable CRNs do indeed have a special structure that allows the relevant Gröbner basis to be computed rapidly, on a time-scale comparable to Gaussian elimination. This is because the fundamental 'building blocks' of the RPA polynomial (i.e. its subsidiary polynomial invariants), always reside in the rowspan of the CRN reaction equations (see Sections S3 and S4), and the RPA polynomial can always be constructed from these by multiplying each by a 'concatenating monomial' (where necessary, depending on the stoichiometric dependence of the invariants, see Section S3), and adding the resulting polynomials. As a consequence, the most computationally-demanding component of the resulting algorithm for an RPAcapable CRN (cf. general polynomial systems) is the computation of these subsidiary polynomial invariants (also known as `complex linear invariants' [8]) from the CRN reaction equations. Note, in particular, that where only ℝ-linear combinations of the rows of ℒ( ) are computed, the S-polynomials of Buchberger's algorithm are identical to the row combinations obtained during Gaussian elimination (applied to the matrix ℒ( )). Gaussian elimination is well-known to be a polynomial-time algorithm, which terminates rapidly even for relatively large numbers of variables (or complexes, in the case of CRNs).
Later, in Section S4, before providing the detailed mathematical development of our approach for computing RPA-relevant polynomials in the rowspan of the rate equations -for low deficiency CRNs ( = 0, 1) in the first instance, and then for CRNs of arbitrary deficiency ( > 1) -we explain that: " … it is only on account of this special 'almost linear' coordinate change that the feasibility of this approach is guaranteed: as we noted in our remarks following Theorem 1, the general problem of computing Gröbner bases, for general collections of polynomials, is known to be NP-Hard. But the fact that the RPA polynomial can always be computed by combining a collection of polynomials in the rowspan of the system means that the algorithm that extracts the RPA polynomial from the CRN reaction equations does not differ substantially from Gaussian elimination. Moreover, even for an arbitrary choice of two projection variables, our statement of Theorem 1 offers opportunities for additional computational efficiency when computing a suitable Gröbner basis insofar as only one particular elimination ideal is required (i.e. that comprising the two chosen variables). Without prior knowledge of this particular property of the RPA polynomial, a full set of all elimination ideals for the system would be required, which necessitates an expensive lexicographic monomial ordering in the execution of the algorithm. With only one elimination ideal required, a more efficient block ordering may be chosen, with a comparatively fast (e.g. degree reverse lexicographic) ordering imposed on all but the two chosen variables." In addition, we provide much more extensive explanations throughout our SI of the mathematical technicalities underpinning this 'almost-linear' coordinate change that can always extract the RPA polynomial from the chemical reaction rates of an RPAcapable CRN. The reviewer is quite right to be astounded by this special property of CRN-capable networks. It is truly remarkable (and entirely non-obvious) that the requisite mathematical structure of chemical reactions that are compatible with RPA should be characterised by such a simplicity as to obviate the computational 'black holes' that can exist for general polynomial systems, including CRN models that do not exhibit RPA.
• Form of the RPA polynomial: The whole approach hinges on the RPA polynomial having the form = ( , )( − ) where ( , ) ≠ 0 and is a non-RPA capable variable which forms the kinetic pair to . It is unclear why ( , ) can only be a function of one additional variable, apart from the output variable. The proof given on page 12 in the Supplement does not satisfactorily explain why the ideal ∩ ℝ ̅ will contain polynomials in and that are not in ∩ .
As an example consider the following network: Furthermore assume that species X2 … X5 participate in several reactions that do not involve X1 but can be catalysed by it. In this scenario the RPA polynomial would be = − . Hence the function depends on two non-RPA variables and . Please explain how this example is consistent with the form of ( , ) stated above.
Our response: Let us first respond to the reviewer's question regarding the proposed network, and its relationship to the two-variable requirement, before returning to the reviewer's concerns about our proof to Theorem 1.
The reviewer proposes two chemical reactions, and adds that we should also 'assume that species X2 … X5 participate in several reactions that do not involve X1 but can be catalysed by it'. Indeed, in the specific reactions proposed by the reviewer, X2 and X3 are only consumed, not produced, whereas X4 and X5 are only produced, not consumed; thus, there must be some additional reactions involving these four variables, so that X2 and X3 can be produced and X4 and X5 consumed, giving the CRN the potential to achieve a steady-state in the positive orthant.
Let us consider several possible ways to 'close' the reaction set (i.e. add in reactions that produce X2 and X3 and consume X4 and X5). First, let us consider two additional closure reactions of the simplest possible form: We emphasize that these are not suggested to be the 'best', or most biologically meaningful, or realistic, closure reactions, only that they constitute a simple solution that allows us to proceed in the first instance. Let , , , , be the concentrations of the five species. In this case, the mass-action equations for this collection of four reactions contain a boundary variable, -a concept which we clearly and prominently define ("Definition 1") at the beginning of Section S1.5 in our SI, in the context of defining what we mean by a variable, before proceeding to invoke this terminology in our statement and proof of Theorem 1 (the two-variable kinetic pairing theorem). Indeed, the rate equations that are induced by this collection of four reactions are: In any case, it is clear from these equations (and, indeed, clear from the corresponding reactions themselves) that is a boundary variable for this CRN, since both and only appear in the monomials of these rate equations in the form of the factor . (In other words, neither species occurs either alone or contributes to a monomial that also contains some other species (e.g. or , say). This CRN thus comprises five species but four variables: { , , , }. But in this particular CRN, the boundary variable technicality happens to be a moot point since this CRN is not capable of RPA (or ACR). This is because < , , , , >∩ ℝ , =< >. In other words, the generator of the ideal consisting of all polynomial consequences of , , , and containing only the variables and is not an RPA polynomial. Thus, for all steady-states of the system, regardless of parameter choices or initial conditions, = 0, and RPA is thereby not imposed on , despite the fact that one of the reaction rates ( in this instance) is an RPA polynomial. In fact, our Theorem 1 makes this important subtlety crystal clear: as counterintuitive as it may seem, it is not enough for an RPA polynomial to be one of the reaction rates, or be otherwise contained in the ideal generated by a CRN's reaction equations. Rather, an RPA polynomial must generate the ideal ∩ ℝ , and this ideal must be principal (i.e. the RPA polynomial must be the sole generator of ∩ ℝ , ) in order for the CRN to exhibit RPA. Just to emphasize, ∩ ℝ , constitutes the projection of the ideal < , … , > onto two variables (again, being mindful of how we define a variable) -the putative RPA variable along with one non-RPA variable.
Although we carefully define and explain the key concept of a boundary variable (see Definition 1 and ensuing discussion in Section S1.5), and elaborate on the concept again in Example 2 in Section S3.2, where we carefully explain that the CRN in question has nine species but eight variables due to the presence of a boundary variable (see discussion after Equation (34)), we recognise that these points must not have been sufficiently clear to the reviewer in our original submission. In our extensively revised SI, we have now provided the reader with a brief overview of the reviewer's CRN example in order to further highlight and explain the concept of a boundary variable in Section S1.5. We have also supplemented our explanations with the alternative terminology "power product" (another word for monomial) just to be absolutely sure that the concept will be absolutely clear to the general reader. We thank the reviewer for providing such a simple and clear example of a CRN that contains a boundary variable to enhance our explanations of this key concept.
Please allow us the opportunity to expand even further on the subtleties that RPA imposes on the detailed underlying structures of CRNs via some additional examples, since the reviewer has been kind enough to offer such a fruitful and illuminating CRN example. In particular, let us now consider two other possible alternatives for closing the reaction set provided by the reviewer. The reviewer specifically suggests a scenario where the added reactions involving X2 … X5 "do not involve X1 but can be catalysed by it". In line with this suggestion, let us suppose that either one, or possibly both, of the two additional reactions include X1 as both a reactant and a product, thereby endowing it with the properties of a catalysing enzyme that is neither produced nor consumed by the reaction. (We assume this is what the reviewer means by not involving X1 "but can be catalysed by it"? The wording is not 100% clear as to what the reviewer actually had in mind here. We apologise if we have misunderstood the reviewer's intentions.) Suppose, then, that the reaction associated with rate constant is now replaced by This altered reaction leaves , and unaffected, but and now become: We observe that is still a boundary variable for this altered CRN. When we project the ideal generated by , , , , (using the updated versions of and ) onto and , we still obtain < >. Thus, the CRN is again thwarted in any attempt to impose RPA on , since = 0 at all steady-states for this CRN (regardless of parameters or initial conditions).
Finally, let us consider the possibility that one of the species that appears in the boundary variable (say , for argument's sake) catalyses one of the additional closure reactions, rather than . In this case, the reaction associated with rate constant is replaced by X2 + X4 → 2X2, and and now become: In this case, is no longer a boundary variable since now also features in the monomial . Thus, this CRN comprises five variables (all species): { , , , , }. So, to test for RPA, we now attempt to project onto x1 and x2. In this case, < , , , , >∩ ℝ , contains only zero (the additive identity for the ring). As expected, the ideal < , , , , >∩ ℝ , also contains only zero. Thus, the CRN remains incapable of RPA. If, instead, we attempt to compute < , , , , >∩ ℝ , , , projecting onto the three variables that appear in the apparent RPA polynomial in the reaction rate, we find that this elimination ideal is given by < >. This confirms that for this CRN, either or (or both) is zero for all possible steady-states of the system. Therefore, this CRN fails to exhibit RPA (at or otherwise).
We warmly invite the reviewer to verify for himself/herself that the three versions of this particular CRN cannot, under any parametric circumstances, exhibit RPA. To assist, we have simulated the second version of the CRN discussed above for 10,000 randomised parameter sets, ∈ (1,100), and initial conditions, (0) ∈ (1,100), see Figure below. From these histograms, the significance of having one rate equation ( ) in the form of an RPA polynomial is clear: If, for a specific random choice of parameters and initial conditions, reaches the value / before or reaches zero, then the -coordinate stops evolving and does attain the 'apparent' setpoint at steady-state. But if either or reaches zero first (before reaches / ), then remains at whatever value it reached at the moment or reaches zero. Thus, the steady-state value of is a matter of chance, dependent upon the vagaries of parameters in comparison with initial conditions. There is a bias in the steady-state value for , with mode of the distribution at / , but can achieve any value in principle at steady-state. Thus, the CRN does not exhibit RPA. The fact remains that our Theorem 1 provides a definitive and universal condition that must be satisfied by any RPA-capable CRN -no exceptions, not even for CRNs that happen to have an RPA polynomial in the rowspan!! (Interesting, the authors of the Cappelletti et al. paper [2] appear to have overlooked this subtlety). We hope that the reviewer can now begin to appreciate that these remarkable and highlyrestrictive conditions on the detailed graph structures of RPA-capable CRNs could never possibly have been gleaned just from our prior work in [1], or by any other prior work (including the work of Cappelletti et al [2]). All three examples we analyse above are incapable of RPA since they violate the conditions specified in our Theorem. Specifically, in each case, the intersection of the ideal generated by the reaction equations with a suitable polynomial ring in two variables (ie. one that comprises one RPA variable and one non-RPA variable, once again being mindful of what we mean by a variable), is not generated by a single polynomial that takes the form of an RPA polynomial. The fact that each CRN induces a reaction rate ( ) that happens to be in the form of an RPA polynomial is (surprisingly, perhaps) irrelevant.
We warmly encourage the reviewer to consider other chemical reactions to add to the pair he/she has proposed. Unlike the Cappelletti et al. [2] study (for example), our study goes far beyond simply 'testing' for the RPA capacity of particular CRN examples (and identifying integral variable(s) if the test is successful). Indeed, our article clarifies in detail (and even more explicitly and incisively now due to our extensive revisions, thanks to the reviewer's helpful feedback and suggestions) how the deficiency of an RPA-capable CRN relates to topological features of the CRN, and concretely delineates the role of this key algebraic invariant in controlling the formation of RPA-relevant 'complex linear invariants' in the CRN's rowspan. For this reason, we can go much further than simply proving that the three specific examples discussed above (encompassing the reviewer's suggested reactions) cannot exhibit RPA: The two reactions proposed by the reviewer cannot impose RPA on no matter which reactions are added to the original two. This is because the two given reactions are linearly independent, and the CRN formed from these reactions has deficiency zero. Moreover, the two complexes that contribute to the apparent RPApolynomial in the reaction rate are not members of the same terminal SCC. It is clear from our discussion of the decomposition of CRNs into algebraically independent subnetworks -a decomposition which is closely tied to the partition of deficiency among independent subsets of reactions (see SI Section S1.3) -as well as our extensive discussions of rowspan polynomials (complex linear invariants) throughout the remainder of our supplement, that no possible addition of extra reactions can turn the original two into an RPA-conferring set. We briefly note in Section S4.4.1 (see also footnote in that section) the special case where a deficiencyzero collection of reactions can engender RPA -namely when the two relevant complexes reside in the same terminal SSC, and correspond topologically to a 'trivial' (isolated) connector node. The pair of reactions suggested by the reviewer are clearly not of this type. No previously published work -not Shinar and Feinberg [7], not Cappelletti et al. [2], nor any other published work on RPA/ACR -is sufficiently general and all-encompassing as to rigorously and definitively delineate the strict requirements for RPA in any CRN -regardless of size or deficiency. We emphasize again, in the most strenuous possible terms, that precisely none of these insights follows straightforwardly from our previous study in [1], which gives absolutely no insight whatsoever into viable chemical reaction structures that could impose RPA on any molecule consumed or produced by those chemical reactions.
Having thoroughly addressed the reviewer's query regarding "how this example is consistent with the form of g(x,y) stated above" we now return to the reviewer's queries about our proof to Theorem 1. The reviewer states that "the proof given on page 12 in the Supplement does not satisfactorily explain why the ideal ∩ ℝ ̅ will contain polynomials in and that are not in ∩ ." Just to clarify, what our proof to Theorem 1 actually says is that, given two variables and that do not exhibit RPA, in addition to a variable that does exhibit RPA, then ∩ ℝ ̅ = ∩ ℝ , , will contain polynomials in and that are not contained in ∩ (the ideal consisting of all polynomial consequences of the rate equations for the CRN that also vanish at = , the setpoint for ). The reviewer is not specific as to why he/she feels that this technicality unclear, but we note first and foremost that should ∩ ℝ , , contain a polynomial in and only (not also including ), then this polynomial is ipso facto not in (since all polynomials in contain by supposition). We also emphasize that the ideal consists of all polynomials in the ring ℝ[ , … , ] that vanish on the variety ( , … , ). In any case, our theorem is fundamentally concerned with which projections (i.e. choices for ̅ ) are guaranteed to be entirely contained within the ideal ∩ .
Our proof to Theorem 1 succinctly points out that if we project the polynomial consequences of the system ( ) onto three variables, two of which ( , ) do not exhibit RPA, then there will be a generator for the associated elimination ideal involving only and (and not -the RPA-capable variable). The reason for the existence of such a generator is that a system whose steady-states can be regulated by any number of independent disturbances (or input stimuli) must be able to adapt to each disturbance individually, corresponding to a single degree of freedom. Once the value of any one non-RPA variable ( , say) is specified (as a result of setting the magnitudes of the various possible disturbances), then the value of any other non-RPA variable (e.g. ) is thereby also specified. The value of , by contrast, is not determined through an alteration to any single degree of freedom (corresponding to one of the independent disturbances) since it exhibits RPA, and is independent of these disturbances. Thus, it follows that some polynomial involving the two non-RPA-variables and is a generator of ∩ ℝ , , .
Now if, for the sake of argument, there were a single generator of ∩ ℝ , , containing all three variables, this would mean that the value of could only be determined once the values of both and were specified. On the other hand, if the generator(s) of ∩ ℝ , , were to contain a polynomial in just , then this would imply that the univariate polynomial either has no positive real roots, or possibly no real roots at all (both of which scenarios corresponding to the non-existence of a steady state for ), or has at least one positive real root whose value is a rational function of the system rate constants (a scenario that corresponds to RPA in ). All of these conditions violate the assumption that is a non-RPA variable -a contradiction.
Of course, Theorem 1 could not possibly be true in general if variables are simply taken to be species. We point out in the preamble to Theorem 1 that in most cases, the variables of a CRN are the species, however should a boundary variable exist, then this must be accounted for in the identification of variables for the CRN (since the component species of the boundary variable cannot be considered variables in that case).
We hope these clarifications have been able to assure the reviewer of the correctness of this Theorem and its proof. We have now added a condensed version of our explanations above as a 'Remark' following the proof to Theorem 1 (and following Definitions 2 and 3). We also decided to remove the reference to ACR in the statement of the theorem (and hence in the proof), and instead added a comment on this point to the list of Remarks following the Theorem. (The reviewer's helpful suggestion to give a clear and explicit definition for RPA at the outset now makes it unnecessary to add in descriptions of ACR, and its relationship to the more general case of RPA, to our statement and proof of Theorem 1. We cannot thank the reviewer enough for such supportive advice). We do, of course, recognise that not all mathematical details of Theorem 1 could possibly be obvious, or clear, to the entire readership of a multidisciplinary journal. It is for this reason that we provide many illustrative examples throughout our SI to demonstrate how the fundamental mathematical principles encapsulated by Theorem 1 are reflected in a range of simple RPA-capable CRNs.
Now to the reviewer's next concern: Moreover, it is not immediately clear how existence of multiple opposer/balancer modules translates into existence of existence of corresponding RPA polynomials. This has not been explained in sufficient detail.

Our response:
We thank the reviewer for this very helpful advice that we make a clearer connection between the existence of multiple topological modules and the corresponding existence of associated RPA polynomials.
Due to our extensive revisions, our manuscript (and especially our SI) now makes this connection very explicit, and explains the underlying technicalities in considerable detail. In particular, our extensive explanations on the decomposition of RPAcapable CRNs into algebraically-independent subnetworks underscores the fact that such independent subnetworks can be analysed separately as to their RPA-capacity (or otherwise), and that the RPA-capacity of 'multi-modular' networks corresponds to the existence of multiple independent RPA-capable subnetworks within such a decomposition. We state explicitly in Section S4.2, for example, as we begin to introduce the mathematical process by which higher deficiency ( > 1) CRNs can be accommodated within our framework, that "… for CRNs for which multiple independent subnetworks contribute independently to RPA, these independent subnetworks correspond to distinct modules from a topological perspective (see [1])." We also emphasize again in Section S4.5, where the general principles are summarised, that "…It is clear from the analysis in the preceding sections that a CRN can be decomposed into lower deficiency independent subnetworks corresponding to the individual modules of a multi-modular network and, in the case of Opposer modules, the controller portion of the module, which can be analysed separately." Our study thus makes it clear that each topological module, which directly corresponds to an algebraically independent subset of the CRN reactions, engenders its own RPA polynomial, specific to that module/subnetwork. Nevertheless, Theorem 1 guarantees that for a multi-modular CRN, choosing one RPA-capable variable (associated with any RPA-conferring module) and one non-RPA-capable variable (from any module, including one that is distinct from the one containing the RPA-capable variable), and projecting the entire system onto those two variables, will yield an RPA polynomial containing those two variables. Thus, Theorem 1 is valid for any choice of two variables with the property that one is RPA-capable and one is non-RPA-capable.
We do acknowledge that in our discussion of the topological consequences of Theorem 1 in Section S3 (see for example Figure S4, and Figure 4 in the main paper), we are referring primarily to RPA in a CRN corresponding to a single module. This is mostly from the point of view of (i) highlighting the existence of an overarching mechanism by which a non-RPA capable variable ultimately imposes RPA on a variable, not necessarily directly; and (ii) distinguishing the (topological) significance of an RPA-capable variable being a regulating, or a regulated, variable.
• Connection to existing works: In [4] the authors consider RPA systems that are maximally robust and find simple linear-algebraic structural conditions that characterise this property in both deterministic and stochastic settings. How do the results in this paper connect to the results in [4]?

Our response:
We are grateful to the reviewer for introducing us to such an interesting and impressive recent article [4], of which we were hitherto unaware.
An obvious point of difference between this article [4] and ours is that these authors also consider RPA in the stochastic setting, while our study is focused squarely on the deterministic setting in the interests of delineating the full space of RPA-capable CRNs in the greatest possible generality. Beyond this, the article by Gupta and Khammash focus on a particular type of RPA they refer to as maximal RPA, or "maxRPA", where 'the setpoint for the key output variable depends on the least number of network parameters, and is insensitive to all the others' -arguing that 'it makes strong evolutionary sense for intracellular RPA networks to attain (such a) maximal robustness'. In this connection, the authors claim to 'work with a stronger notion of robustness' and that their 'characterising conditions for RPA networks are more restrictive than the topological requirements expounded in (our previous work [1], for example)'. The study is able to identify, concretely, a number of mathematical conditions required for the existence of the maxRPA property, which includes (among other requirements) that the setpoint should be a function of (only) two biochemical parameters. In any case, the authors succeed in their goal to 'mathematically characterise and study the structure of biomolecular networks that constitute such maximal RPA (or maxRPA) systems' -a most impressive achievement -and propose a novel internal model principle (IMP) pertaining to these maxRPA networks.
This focus on maxRPA networks alone constitutes the major distinction between [4] and our study, since our goal is to characterise all possible RPA-capable CRNs -as opposed to specific classes of RPA-capable CRNs (e.g. the maxRPA class [4], the deficiency-one Shinar-Feinberg class [7], the class containing a 'linear constrained integrator' [2], etc). As we noted earlier, this generality is important in the context of understanding, at an abstract level, the ability of RPA-capable networks to cope with new disturbances (i.e. those that were not present in the original CRN -e.g. the addition of small-molecule enzyme-inhibitors). This is particularly crucial for signalling networks for which the detailed structures of the underlying CRNs are unknown (in the context of signal transduction in metazoan cells, for instance). Moreover, we respectfully question the authors' point of view on the suggested evolutionary advantages of a setpoint involving only two network parameters. We contend that it is not the number of parameters defining the setpoint, but the nature of the parameters that is most significant to the functional robustness of the system, and ultimately, the survival of the cell/organism. A CRN could have a setpoint determined by a hundred parameters, but if all these parameters are 'stable' in the sense that they are determined entirely by the chemical properties of the interacting molecules (e.g. association/dissociation constants, or catalytic constants), alterable only by mutation or by factors such as temperature that vary very slowly on the timescale of adaptation, then this may still constitute an entirely useful form of RPA from a biological standpoint. By contrast, a CRN could have a setpoint determined by only two parameters (the minimum, as in the antithetic integral control motif), but if these are highly labile parameters relating to, for example, the rate of synthesis of a particular molecule, or the total expression level of a particular molecule, then this could potentially be problematic.
Consider, for instance, the deficiency-two Shinar-Feinberg model of the EnvZ-OmpR motif. The setpoint for pOmpR due to this CRN involves eleven parameters (calculated from two independent setpoints arising from the balancer and connector contributions). But all of these parameters are either association/dissociation constants for pairs of molecules, or catalytic constants, and are thus highly stable in magnitude. By contrast, for the antithetic control mechanism, the setpoint for the sensor molecule involves just two parameters. But both of these parameters capture the rate of synthesis of an antithetic molecule which could, in principle, vary quite significantly on the timescale of adaptation (depending on the regulatory mechanisms governing the transcription and translation of the molecules in question).
Moreover, perturbations to the setpoint that are either stable (e.g. due to a mutation in an RPA-relevant signalling protein) or vary on a timescale that is significantly longer than the timescale of adaptation, may not alter the functionality of the biological system whatsoever. One of the most widely-discussed examples of unaltered biological function in the face of a stably varying setpoint is the robustness of tumbling frequency, and precision of adaptation, in E.coli chemotaxis (see Alon et al. [12]). In that study [12], genomic alterations that affected the expression levels of various chemotaxis proteins were introduced, which stably altered the setpoint for the steady-state tumbling frequency as well as the adaptation time. But exact adaptation was still observed (albeit with an altered setpoint), and the chemotaxis response of the cells was unaffected. In this context we recognise that, for singlecelled organisms, significant temperature variations could certainly alter biochemical parameters; but provided as temperatures vary slowly in comparison with adaptation (RPA), this might not affect the functionality of the organism at all. Moreover, aside from what type of setpoints are most evolutionarily advantageous, we respectfully point out that the RPA responses frequently observed in highly complex signalling networks, such as cancer signalling networks in human cells (on the basis of time-course data -see [13], for instance) are not necessarily attributable to maxRPA networks; there is simply no evidence to support this. In fact, it is currently entirely unclear what CRN architectures are generally responsible for RPA in complex signal transduction networks (in either normal or diseased states).
We hope the reviewer will not object to a friendly and respectful debate on these points, in the interests of moving scientific frontiers forward. We contend that, to make sense of the highly complex signalling networks that arise in nature, it is absolutely essential to have access to completely general description of the mechanisms by which RPA could be implemented in CRNs. This is the brass-ring objective that we achieve in the present work.
In raising all these points, we are by no means dismissing the study by Gupta and Khammash [4], which is incredibly interesting, and extremely impressive (mathematically and otherwise); we only point out that it is not completely general vis-à-vis RPA-promoting CRN structures (and makes no claims to be, since its focus is on the particular class of RPA-capable CRNs known as 'maxRPA'). It is a beautiful paper, and we are proud to include it in our reference list in the main manuscript.

Concluding Comments to Reviewer #1:
We sincerely thank the reviewer again for such a generous investment of time and painstaking effort, and for his/her many supportive and helpful suggestions for making our scientific goals, and our findings, much clearer and more explicit. We are truly grateful to have this opportunity to defend our work, and to provide an extensively revised version of our manuscript that offers a much more comprehensive exposition of all the relevant supporting technicalities, and a much more compelling presentation for a general scientific audience.

Response to Reviewer #2
Araujo and Liotta develop a set of precise criteria that a system needs to fulfil in order to perfect adaptation. This is remarkable and an extremely valuable contribution to the literature in three ways: (i) the search for "design principles" is taking centre stage in synthetic (but also e.g. developmental systems) biology; the first author has proven a leader in distilling precise mathematical criteria for robust perfect adaptation, and I hope that this study will inspire more work in this area, especially in considering the mathematical properties of chemical reaction networks to guarantee different types of behaviour; (ii) in Mathematical biology there are few instances where the mathematical statements can be so precise that biology will have to "obey" these statements. This is one such instance. Finally (iii) the maintenance and control of robust adaptation is of considerable and far reaching biological relevance. I like this manuscript for each of these three points a lot.
The most interesting mathematical aspects of the work have, unfortunately but predictably been relegated to the supplementary information. The results in Figure 2, for example, are very clearly explained and easy to follow in the SI. I personally enjoyed the SI a lot, and while the discussion e.g. in lines 207-215 is clear, it may be too terse for some readers to follow. The SI by contrast was very clear. I think it would help readers to understand how translatable to other phenotypes this type of analysis is, or if RPA is particular in this regard of allowing such general mathematical statements to be derived.

Our response:
We are truly indebted to this reviewer for his/her incredible generosity in considering our work in such careful detail, and for such genuinely helpful feedback and thoughtful suggestions for further improvement. We are immensely grateful to have the benefit of an 'extra set of eyes' cast over our work, to identify additional opportunities to make our exposition more compelling. We have noted all revisions to the main manuscript, and to our SI, with red text, to make it easy to identify which parts of our article have been updated since the original submission.
Regarding the issue of relegating much of our detailed mathematical development to the Supplementary Information, it is true that the severe word limits in Nature Communications (6000 words for main text, excluding figure captions) have forced us to limit our exposition in the paper itself to an accessible overview of the major conceptual innovations of the study, with extensive references to all the technical details in our Supplementary Information. A significant additional challenge stems from the fact that solving the RPA problem in complete generality, at the level of intermolecular interactions, has required us to draw together a number of distinct mathematical languages -chemical reaction network theory (CRNT), engineering control theory, and algebraic geometry -into a unified framework.
We greatly value the reviewer's insightful suggestion that "it would help readers to understand how translatable to other phenotypes this type of analysis is, or if RPA is particular in this regard of allowing such general mathematical statements to be derived." We now include following text at the conclusion of our article: The quest to uncover the fundamental 'design principles' that constrain complex signalling networks in nature to implement biologically important functions is considered to be one of the most important and far-reaching 'grand challenges' in the life sciences [1][2][3][4][5][6]. On the basis of the present study, along with our earlier topological study at the network macroscale [7], RPA currently stands alone as a keystone biological response for which there now exists a universal explanatory framework -one that imposes strict and inviolable design criteria on arbitrarily large and complex networks, and one that now accounts for the subtleties of intricate intermolecular interactions at the network microscale. These universal RPApermissive design principles now represent a launching-point for future explorations of more complex phenotypes -including some classes of embryonic patterning problems, for instance, where integral control is known to play a role in promoting adaptation of segmentation boundaries to variations in organism size [8,9]. The identification of universal design principles for many other complex phenotypes, such as Turing patterns [4,10,11] and multistability/switching-responses [12,13], is likely to prove more challenging due in part to the central role of equilibrium stabilities, or instabilities, in these responses. These 'grand challenges' remain open, and we hope that our study will inspire bold new mathematical thinking in these vitally important directions.
Small points: -much of chemical reaction network theory is obscure to most readers. I would like to see a clear and easy definition of integral control that is accessible to non-expert audiences. Deficiency is maybe another such concept and it would be good to define it on line 146.

Our response:
We truly appreciate such helpful suggestions to improve the clarity of our exposition. In our newly revised manuscript, we give a thorough yet accessible explanation of the concept of integral control, and the internal model principle (IMP), via the introduction of a new Figure in our Introduction section ( Fig. 1 in our revised submission) -see below. We have also carefully edited our explanations of integral control in-text to complement this new Figure, and to ensure that the concept (and its relationship to the goals of our study) are clear to the readership of a multidisciplinary journal.
We have also vastly expanded our clarifications on the concept of deficiency, especially in our Supplementary Information where we give far greater detail on the concept where we first define it in SI Section S1.2, and particularly its relationship to the flow of biochemical information (and ultimately topology) in RPA-capable CRNs. We have also added in an entirely new section (SI Section S1.3) on the key concept of algebraically-independent subnetworks of a CRN (which we didn't discuss until much later in the Supplement, in our original submission), and clearly outline the role of deficiency in partitioning the reactions of a CRN into such independent subsets. In the main paper, we also introduce the concept of deficiency much earlier than previously, and point out early that this is a key algebraic invariant which has fundamental (and previously unrecognised) consequences for the implementation of integral control by CRNs. Later in the main paper, where we give further detail on the underlying mathematical technicalities, we give a clear definition of deficiency, and give an example of how to calculate it, in Figure 5.

b) A suitable coordinate change should be able to recast an RPA-capable system into integral feedback form, even if there is no feedback present in the network. As shown, a linear transformation is sufficient to identify an output-driven internal model for this particularly simple incoherent feedforward motif (Balancer module); = / (setpoint) at steady-state. (c) A model that employs feedback is frequently simpler to recast in 'integral feedback' form, with an output-driven internal model; here = (setpoint) at steady-state. Note that the reaction rates selected for illustrative purposes in (b) or (c) cannot be induced, under the law of mass action, by any CRN 20 .
-The authors show that their criteria hold for all instances of RPA. I was wondering how easily in practice these mechanisms could be lost? Is it easy to identify points that disrupt RPA?
Our response: With access to a universal and complete solution to the RPA problem in CRNs, we can now delineate, precisely and definitively, the circumstances under which various classes of disturbances, or network alterations, can either preserve the RPA property or cause RPA to fail. We have not commented on this matter in detail in our manuscript, as this is the subject of ongoing work. One of the most interesting classes of disturbances against RPA-capable CRNs is the addition of enzyme-inhibitors of various types (competitive, non-competitive, mixed, uncompetitive, etc.). Addition of one more enzyme-inhibitors to a CRN, whether highly specific to its target or not, corresponds to the addition of extra linkage class(es) to the CRN. We now know that certain inhibitor mechanisms are deficiency preserving, while others aren't, and their effects on the RPA property (depending on the relationship of their target protein(s) to the overarching structure of the CRN) are complex and counterintuitive, and in some cases very surprising. Generally, the RPA property is very difficult to perturb via molecular-targeted inhibitors, although we have been able to identify very specific and definitive conditions in which this is possible. Other classes of disturbances (such as those we itemize in our newly revised Supplementary Information, at the beginning of Section S1) are much easier to study from the point of view of altering the RPA property. In any case, it is only by having a completely general mathematical description of RPA-permissive reaction structures, that accounts for all possible RPA-capable CRNs, that we can make such definitive conclusions on the preservation or destruction of the RPA property in the abstract. This is particularly valuable in the context of particularly complex signalling networks (e.g. in cancer signal transduction) for which detailed CRN structures are unknown (and possibly will never be known).
In any case, we are planning to publish this work (separately) very soon, and look forward to sharing this.
-there is a vast literature on design principles which could be touched upon at least in passing (limit cycles, multistability, switch-like behaviour, Turing patterns), especially if there is scope for applying similar concepts in these contexts.
Our response: This is such a wonderful suggestion, which we greatly value. As we noted in our first point above, we have added an extra paragraph that comments on this very issue at the conclusion to our manuscript. "The quest to uncover the fundamental 'design principles' that constrain complex signalling networks in nature to implement biologically important functions is considered to be one of the most important and far-reaching 'grand challenges' in the life sciences [1-6] …. Etc."

Concluding Comments to Reviewer #2:
We cannot thank the Reviewer enough for such a generous investment of time and effort, for such incredibly supportive feedback, and for so many insightful suggestions for additions and clarifications to our exposition. We are truly grateful.

Robyn P. Araujo and Lance A. Liotta
Overview: The aim of this paper is to provide an algebraic characterisation of the hidden integral controller and construct an algebraic procedure to identify it in adaptationcapable networks. It is shown that this procedure is intimately connected to the structure of the networks via the well-known deficiency theory for chemical reaction networks.
Recommendation: I am grateful to the authors for revising the manuscript to address many of the points that I raised in the first review. Even though the paper has significantly improved, many of the more serious concerns remain. These concerns, which are outlined below, should be fully addressed before the paper can be reconsidered again for Nature Communications.
1. Connection to integral control: The paper claims to uncover universal structures for embedded "integral control". However in order to successfully demonstrate this, a coordinate transformation (possibly nonlinear) needs to be constructed that gives rise to the "integrator" variable. It does not matter if this coordinate transformation is done in a single step or in multiple steps (as proposed by this paper), but its existence needs to be shown in order to identify the integral mechanism. In the Supplement the authors write We note that, although there should always exist some single nonlinear coordinate change that yields the RPA polynomial in principle, for any mass action system, identifying such a single coordinate change in practice may be extremely difficult in all but the very simplest RPA-capable CRNs.
I agree that proving the existence of an integrator may be difficult, but unfortunately it is unavoidable for showing integral control. Let us look into this issue in more details.
Suppose the dynamics is governed by a system of ODEs: where x is the n-dimensional state of the system and f (x) = (f 1 (x), . . . , f n (x)) are the rates of evolution of all the species. In order to exhibit an integrator, we need to find a real-valued function H(x) such the dynamics of z(t) = H(x(t)) is given by the ODE: where ∇ is the gradient operator, • is the usual dot product and ρ(x) = p(x)(x i − c) is what the authors call the "RPA polynomial" for the RPA variable x i . Suppose one is able to find nonlinear polynomials h 1 , . . . , h n (as in Theorem 1 of the paper) such that However this does not necessarily mean that a function H(x) can be constructed so that (0.1) holds. Observe that such a function H(x) would need to have the following partial derivatives Unless each h i (x) is a constant (in which case H(x) is linear), it is unclear why such a function H(x) would exist. Since the approach in the paper uses nonlinear elimination steps (via concatenating monomials), the existence of H(x) cannot be ascertained. In fact if for a pair j, k we have then certainly H(x) cannot be constructed with partial derivatives h 1 , . . . , h n . This discussion shows that while existence of a RPA polynomial may be necessary for integral control, it is not sufficient and hence this property does not characterise integral controllers (a main claim of the paper). This gap between necessity and sufficiency can be bridged if the authors can show that for RPA networks one always has 2. Questions about the Kinetic Pairing result: In the previous review round, I had proposed a possible RPA network and asked the authors how it satisfies the Kinetic Pairing result (Theorem 1). The network has two reactions given by and there may be several other reactions involving species X 2 , . . . , X 5 that can be catalysed by X 1 but they do not change X 1 . In the rebuttal letter, the authors consider certain instances of such a network, and claim that in these networks all possible steady states lie on the boundary on the positive orthant (e.g. x 2 or x 3 is zero) and hence the networks are not RPA. In particular, the authors state in their rebuttal letter that: The two reactions proposed by the reviewer cannot impose RPA on X 1 no matter which reactions are added to the original two. This is because the two given reactions are linearly independent, and the CRN formed from these reactions has deficiency zero. Moreover, the two complexes that contribute to the apparent RPA-polynomial in the reaction rate f 1 are not members of the same terminal SCC. It is clear from our discussion of the decomposition of CRNs into algebraically independent subnetworks a decomposition which is closely tied to the partition of deficiency among independent subsets of reactions (see SI Section S1.3) -as well as our extensive discussions of rowspan polynomials (complex linear invariants) throughout the remainder of our supplement, that no possible addition of extra reactions can turn the original two into an RPA-conferring set. We briefly note in Section S4.4.1 (see also footnote in that section) the special case where a deficiency zero collection of reactions can engender RPA namely when the two relevant complexes reside in the same terminal SSC, and correspond topologically to a trivial (isolated) connector node. The pair of reactions suggested by the reviewer are clearly not of this type.
To test the authors' claim above, I tried to come up with instances of such networks (i.e. (0.2) + some reactions) which are RPA and yet do not have this issue of steadystate not being in the positive orthant. Please explain how the Kinetic Pairing result fits these networks.
I start with a simpler version of this network with only three species (i.e. species X 4 and X 5 are absent) where I add inflow and outflow for species X 2 and X 3 . So the overall network becomes This network has no boundary variables and so the variables are the same as species. It seems that this network will indeed show RPA for x 1 with set-point c 1 /c 2 . In Figure 1 I plot the simulated dynamics (rescaled to have the set-point of x 1 as 1) for three randomly chosen values of initial states and rate constants, and one can see that RPA holds.
Next I modify this network to have the production of X 2 and X 3 catalysed by X 1 . So the new network becomes This network also appears to be RPA as shown by the simulations in Figure 2.
Finally, I consider another variant of this network where we have species X 4 and X 5 that reversibly bind to each other to produce an inactive compound, and they catalyse the production of X 2 and X 3 respectively. Therefore the new network becomes This network also seems to be RPA, as shown in the simulations in Figure 3.
These networks do not appear to satisfy the RPA characterisation given by this paper. Also consider Example 6.5.3 in [5] which does not seem to fit this result either. These counter-examples cast doubt on the correctness of the Kinetic Pairing Theorem on which the entire paper rests.
3. The set-point may not be a rational function of parameters: In many places in the main paper and the supplement, the system's set-point is said to be a rational function of biochemical parameters. While this holds for most examples, consider the following birth death network The dynamics is given by and so the set-point is not a rational function of c 1 and c 2 ? For more examples, see Examples 2.11 and 2.12 in [3]. Please revise as necessary. The Gröbner basis algorithm to find the RPA polynomial may not terminate. It is mentioned that failure to terminate for a chemical reaction network (CRN) is a prima-facie evidence that the CRN does not exhibit RPA. However this is not mathematically shown. In any case, checking for nontermination of a method is impractical.
In response the authors state in their rebuttal letter that the termination is guaranteed because "all ideals of a polynomial ring (or any Noetherian ring, for that matter) are finitely generated a result of central importance in algebraic geometry, formalized by the Hilbert Basis Theorem". While I agree that termination is guaranteed, I must point out that the reason for my confusion and for raising this termination issue is the following excerpt from the previous version of the Supplement (submitted in round 1): Failure of the Buchberger (or other Grobner basis-computing) algorithm to terminate could thus be adduced as prima facie evidence that the CRN under consideration does not, in fact, have the capacity for RPA.
It seems that in this statement that authors meant "terminate in practical time". I understand that due to the "almost linear" nature of coordinate change, this algorithm would work well for RPA networks. However, does the final RPA polynomial produced by the method depend on the monomial ordering or other choices made by the method? More importantly, if a network is not RPA, how many steps would the method need to confirm this non-RPA property? Such practical considerations must be discussed in the main text, and they are crucial for applying these ideas for characterising RPA in high-dimensional networks.

5.
Connection of examples to existing works: The authors start the Results section with "two simple examples that have eluded all previous systematic methods to detect RPA". Please specify which systematic methods are being referred to here. Also, in the example on Figure 2 it should be mentioned that X 3 is maxRPA and it can be checked from the characterisation result in [2]. Secondly since there is the reaction which does not involve X 1 I do not understand why the term k 6 X 2 X 3 does not enter the expression for d(X 1 − X 2 ) dt in Figure 2C. Please check. If the term k 6 X 2 X 3 is present then it cannot be eliminated with the concatenating monomial O 1 . On the other hand, if this reaction is removed (i.e. this term is absent) then the overall network simply becomes a trivial RPA network where the output of one RPA network (i.e. the network comprising X 3 − O 1 ) is passed as an input to another RPA network (i.e. the antithetic network with X 1 − X 2 ). It is straightforward that connecting RPA networks in series (with catalytic reactions) would still result in a RPA network. Such examples are not appropriate for demonstrating the novel results in this paper.
The second example in Figure 3 seems to be taken straight from [4] (see Fig. 2). This should be clearly stated when the example is being introduced in the main text and also in the caption of Figure 3. Also mention that the linear invariants shown in Fig.  3 can be deduced from the approach in [4] (this is stated in passing in the conclusion but it should be stated more prominently when the example is being discussed.) 6. Many claims without proofs in the Supplement: The Supplement has been considerably revised, but still many arguments are unclear because proper proofs have not been provided or referenced: • Why should eq. (8) hold when eq. (9) holds? Can the networks always be partitioned this way? Please explain.
• On page 21 in the Supplement it states that: Since a perturbation to the CRN that alters the steady state of x j will also alter the steady-states of other non-RPA capable variables (eg. x m ), If R[x] will contain polynomials in x j and x m , and that are not contained in I f ∩ I p .
Why should such a perturbation always exist?
• Also on page 21 in the Supplement it says that The set x now contains two independent (uncoupled) variables in the sense that a perturbation to the CRN that alters the steady-state of one of the variables does not affect the steady-state of the other.
Why does this hold? Please elaborate.
• In general, in the proof of the Kinetic Pairing Theorem the authors work over the ring of polynomials over species-variables x 1 , . . . , x n . Shouldn't the system parameters (i.e. rate constants) be included in this ring, as they would appear in the RPA polynomial? This inclusion of parameters is there in the Singular code but not in the proof. However simply adding the parameters in the ring is not sufficient as the n-th roots of the parameter would be added, as the examples mentioned in point 3 show.

Other minor issues:
• The definition of RPA must be shifted to the main text due to its centrality in understanding the message of the paper.
• Why is the variable x i missing in g(x j ) in figure 4 (main text) and figure S4 (supplement)?
• Replace "consistutes" with "constitutes" on line 72 in the main text. Please run a spell check.
• The paper says that if a network is RPA the integrator is guaranteed to exist. For example the following on page 22 in the main text In principle, there should always exist some single nonlinear coordinate change to extract a single output-driven internal model (Fig. 9a) from systems rate equations, corresponding to a single integral of the systems tracking error (Fig. 9b) or the following on page 16 in the Supplement All classes of RPA, including ACR, thus require some form of integral control.
Please explain which version of IMP can be used to verify this existence. See [1] for a recent review on IMP.
• On page 20 the authors state that "It is striking to note that the original form of the CRN (Fig. 8a) eludes the Shinar-Feinberg theorem, even though the CRN exhibits ACR and has a deficiency of one." However does it fit the results in [4]? Please comment on this. Also explain why the authors found the method for finding linear invariants "ad hoc" (lines 568-569 on page 25).

• On the Supplement page 15 it is stated that
Mass-conservative CRNs therefore have no external stimuli or inputs, and can only be perturbed by altering the total abundances (or concentrations) of the constituent molecules -i.e. by altering the initial conditions.
Why cannot the perturbation come in the form of parameter variation, e.g. of a conversion reaction (that conserves mass).

Reviewers' Comments:
Reviewer #1: Remarks to the Author: See attached.
Reviewer #2: Remarks to the Author: The reviewers have comprehensively addressed my questions and points. This paper addresses a fundamentally important question in biology and it does so using a powerful mathematical framework that the authors have developed and fine-tuned.
The particular appeal of this approach is that it sets a template (and I imagine that in the future we will see many more examples of this type) is that the mathematical analysis has a level of generality that means that the authors results are robust and reliable. Having solid mathematical foundations for a complex biological phenotype is still rare, but the current paper shows that the search for such theoretical underpinnings can be fruitful.
Overview: The aim of this paper is to provide an algebraic characterisation of the hidden integral controller and construct an algebraic procedure to identify it in adaptation-capable networks. It is shown that this procedure is intimately connected to the structure of the networks via the well-known deficiency theory for chemical reaction networks.
Recommendation: I am grateful to the authors for revising the manuscript to address many of the points that I raised in the first review. Even though the paper has significantly improved, many of the more serious concerns remain. These concerns, which are outlined below, should be fully addressed before the paper can be reconsidered again for Nature Communications.
Our response: Once again, we are truly grateful for the huge investment of time the reviewer has clearly made in considering our work in such careful detail. We have reflected thoughtfully on all the reviewer's points and respond to each in turn below.
All changes made to our paper and Supplementary Information (SI) are noted in red text.

Connection to integral control:
The paper claims to uncover universal structures for embedded "integral control". However in order to successfully demonstrate this, a coordinate transformation (possibly nonlinear) needs to be constructed that gives rise to the "integrator" variable. It does not matter if this coordinate transformation is done in a single step or in multiple steps (as proposed by this paper), but its existence needs to be shown in order to identify the integral mechanism. In the Supplement the authors write We note that, although there should always exist some single nonlinear coordinate change that yields the RPA polynomial in principle, for any mass action system, identifying such a single coordinate change in practice may be extremely difficult in all but the very simplest RPA-capable CRNs.
I agree that proving the existence of an integrator may be difficult, but unfortunately it is unavoidable for showing integral control.

Our response:
What we show in this paper is that a particular type of nonlinear transformation is always able to decompose an RPA-capable set of CRN reactions into an organized collection (a 'topological hierarchy') of linear coordinate transformations. Each of these linear coordinate transformations corresponds to a subsidiary integral controller within the CRN. This is a fundamentally different type of integral control implementation from the 'conventional' interpretation of integral control in engineering control theory. Indeed, the nonlinear transformation we propose is not a coordinate transformation. Rather, it is a map from the system reaction rates to a function space with the same algebraic structure as the CRN reactions (i.e. that of a commutative ring), which offers a universal description of the coordinated interactions of (integral computing) internal models -obtained via linear coordinate changes -that holds for all possible RPA-capable CRNs.
We do acknowledge that in our previous rebuttal letter (as opposed to our paper), we had referred to the mathematical transformation in several places as a coordinate change, which was not correct. We apologise for unwittingly sidetracking the discussion through such sloppy (and incorrect) wording in our previous response. We have now made absolutely certain that this incorrect terminology does not appear anywhere in our paper or in our accompanying SI. The reviewer is entirely correct to point out that the nonlinear transformation we present in this paper is not a nonlinear coordinate change.
The conventional view of the 'internal model principle' (IMP) in engineering control theory (particularly for application to problems in systems biology and bioengineering [1,8,9]) suggests that, for robust rejection of constant disturbances, there must exist some decomposition of the system (obtained via a coordinate transformation, if needed) into two distinct subsystems: (i) an 'output-driven internal model' (of the constant disturbance), and (ii) the remainder of the system ( [8,9]). This decomposition will generally require a nonlinear coordinate transformation, although many simple CRNs have been identified (as we point out in our paper and accompanying Supplementary Information) for which a linear coordinate transformation is sufficient. In particular, if the system at hand is not already in feedback form, then the internal model-identifying coordinate transformation must be able to recast the system into an integral feedback form (see references [1] and [9] below, as well as Figure 1b in our paper, for some simple examples of this requirement). If a network can indeed exhibit RPA, then according to this version of the IMP (see [8] and [9], for example, as well as Section 3 in [1]), there should always exist some (single) coordinate transformation that recasts the system into such a feedback system, with an embedded output-driven internal model -even if the system is topologically feedforward structured. This is the concept we are referring to in the excerpt from our Supplementary Information quoted by the reviewer above ('… there should always exist some single nonlinear coordinate change that yields the RPA polynomial … for any mass action system').
Note, incidentally, that (in contrast to the formulation we present in our paper) this conventional interpretation of the IMP ( [8,9]) requires that certain properties of the disturbances be known a priori, including which model variables can be directly regulated by the disturbance, and also requires the prior identification of an 'output' for the model (which must be able to 'adapt' to, or reject, those disturbances). For an RPA-capable system, being able to identify such a nonlinear coordinate transformation, if it exists, simply confirms to us that the system can indeed exhibit robust asymptotic tracking of the output variable's 'setpoint' when subjected to the specified (constant in time) disturbances. It also identifies the specific transformed variable that actually computes the integral of the 'error' (where the error in question is the difference between the instantaneous value of the output variable and its setpoint). But if we already know (using some other analytical method, say) that RPA obtains for certain variable(s) in a specific CRN, then explicitly identifying the nonlinear coordinate change that transforms that particular system into integral feedback form does not actually provide any further useful information. The fundamental point we make in our Discussion is that there is no general way to propose a single nonlinear coordinate change that can reveal the fundamental mechanisms by which all RPA-capable CRNs actually implement RPA, and thereby provide a concrete characterization of the entire solution space to the RPA problem for CRNs.
Consider, for example, Shinar and Feinberg's deficiency-two model of the EnvZ-OmpR osmoregulation network (first presented in the Supplementary Materials of [10], and considered in Figure 3 of our main paper). For this specific network, the existence of RPA (and ACR more specifically) in the molecule phospho-OmpR can be shown using a variety of methods. Shinar and Feinberg, in their Supplementary Materials [10], show through a sequence of manual algebraic substitutions that phospho-OmpR has a steady-state value that is a rational function of the CRN parameters, independent of the total abundances of the various interacting molecules. Karp et al. [6], on the other hand, demonstrate how 'complex-linear' invariants may be computed, and show that for this particular CRN, two linearlyindependent such invariants may be combined to demonstrate RPA in phospho-OmpR, with the same setpoint as previously calculated by Shinar and Feinberg [10]. Now, since this CRN exhibits RPA (at phospho-OmpR), undoubtedly, there must undoubtedly exist some nonlinear coordinate transformation that could recast this system (a 'balancer' module, as we show, with a feedforward structure) into an integral feedback system with a single output-driven (i.e. phospho-OmpR-driven) internal model, where a single (transformed) variable computes the error in phospho-OmpR in comparison with its (known) setpoint. But even if we were to identify the requisite coordinate transformation, explicitly and analytically, for this particular CRN, what would this coordinate map actually tell us (vis-à-vis the RPAcapacity of the CRN) that we didn't already know? More importantly, what could this specific coordinate map possibly tell us about the space of all possible RPAcapable CRNs, or the general (universal) properties of collections of chemical reactions that can confer RPA on a subset of the interacting molecules? Absolutely nothing! The integral-feedback-recasting nonlinear coordinate transformation for an RPA-capable CRN is unique to that particular CRN. In other words, the single nonlinear coordinate transformation required to identify an internal model in the deficiency-two Shinar-Feinberg EnvZ-OmpR model discussed above, if such can be explicitly identified, will necessarily be a different nonlinear map from the one required for the Cappelletti et al. [11] toy model (Example 1 in Section S3.1 of our SI) -despite the fact that the two models are fundamentally alike: both are Balancer modules (topologically speaking) for which a single balancer invariant and a single connector invariant are obtained via linear coordinate transformations, and require a single concatenating monomial to obtain an RPA polynomial.
By contrast, our nonlinear transformation (Eq. 6) is of a fundamentally different character. Nowhere in our paper do we claim that this nonlinear transformation is a coordinate transformation. Importantly, this nonlinear map relaxes the feedback requirement of the IMP in its conventional interpretation (as discussed above, see also Section 6.2 in [9]), and instead preserves the underlying topological structure of the network. Thus, RPA-capable CRNs of Balancer type retain their feedforward characteristics under this transformation. We maintain that referring to an embedded integral control is entirely justified since the nonlinear map in question is always able identify, for any RPA-capable CRN, a topologically organized collection of linear coordinate changes, each one of which identifies a subsidiary internal model which does recapitulate the dynamical structure of the disturbance (i.e. constant-in-time), and which thereby imposes RPA on some characteristic of the network, with its own setpoint. Each such linear coordinate change thereby corresponds to subsidiary control problem with its own linear integral variable.
For the Cappelletti et al. [11] toy model, for instance, the nonlinear transformation which projects the system onto two variables (A and B -see Singular Code in S5.3 of our Supplementary Information) can be decomposed into two linear coordinate changes, corresponding to two key invariants. One of these is a 'balancer' invariant, Together, these two independent computations, each with a linear integral variable, and each conferring RPA on some (topologically important) feature of the CRN, confer RPA on the molecule . The setpoint for is a combination of the setpoints from the two contributing subsidiary (linear) problems. Of course, many CRNs will have much more complicated invariants than these, involving more than two terms per invariant, but the same fundamental principles hold. (This is the universality of the framework we present).
The feedforward model presented by Bin et al. [1] provides a nice illustration of the distinction between the conventional interpretation of the IMP and the one we develop here, as applicable to RPA in CRNs. In Section 3 of [1] the authors review a very simple two-variable incoherent feedforward loop (IFFL) model (Eq. 25) that has been widely studied in the systems biology literature (including in [9]). The two model equations are The authors show that a nonlinear coordinate map that transforms these two equations into a partitioned form, corresponding to an integral feedback structure, is This transformed system has an output = & , and has the desired internal model form, where ! is the transformed variable that integrates the error (in & = ). This nonlinear coordinate change nicely confirms that RPA is possible at , and that its setpoint is / . But a different nonlinear transformation of the system (not a coordinate transformation) is also possible: Unlike the coordinate transformation considered previously, this particular transformation underscores the fact that this simple two-variable reaction system is fundamentally of the same type as both the Cappelletti et al. [11] toy model (Example 1 in our Supplementary Information S3.1) and the deficiency-two Shinar-Feinberg EnvZ-OmpR model we discussed above: a Balancer module, with two independent (linear) invariants -a balancer invariant, and a connector invariantcombined via a single concatenating monomial (in this case, ). This transformation also reveals the setpoint, = / . Of course, we do acknowledge that the rate equations above are not polynomials, the only consequence of which, in the context of our paper, is that we cannot use the powerful automated algorithms of algebraic geometry that compute Gröbner bases (which pertain to polynomials); nevertheless, rational functions are & -smooth functions on (0, ∞), and thus possess the algebraic properties of a ring as required for the application of our Theorem 1. The RPA polynomial so obtained involves two variables (one RPA-capable (i.e. ) and one non-RPA-capable (i.e. )) and is thus the single generator of the relevant ideal. By contrast, the nonlinear coordinate transformation discussed in [1], as noted above, is unique to this particular set of rate equations, and provides no insight into the general properties of all RPA-capable CRNs (or even its fundamental connection to the RPA-conferring properties of the Cappelletti et al. [11] or Shinar-Feinberg [10] models discussed above). Shoval et al. [9] acknowledge that the coordinate transformations that recast systems into integral feedback form 'may well be merely a mathematical construct with no biological meaning' (see Section 6.3 in [9]).
Again, we maintain emphatically that this is a fundamentally different way of looking at integral control from the IMP as it is usually understood in engineering control theory. This approach identifies where, within the collection of reactions themselves, the computations relevant to RPA actually occur at the level of the interacting molecules of the CRN. Identifying a (single) nonlinear coordinate change that transforms the system into an integral feedback system, with a single embedded internal model, can only confirm that RPA will obtain; it does not explain how it is implemented by the original variables of the system (ie. the molecular concentrations).
Being able to characterize the entire space of all possible RPA-capable CRNs in such a general way has enormous practical implications because it provides a completely comprehensive and universal understanding of how to either destroy or maintain the presence of RPA -through evolution, or via experimental or clinical interventions, where new mutational events or exogenous enzyme inhibitors exert their effects at the molecular level. In particular, we can destroy the RPA property in a CRN by disrupting any one linear integral controller. By contrast, retaining the RPA property requires that all linear integral controllers already present be preserved. Cappelletti et al. [11], for instance, consider how a linear integral controller can be preserved when adding in a new (exogenous) stimulus to the network (see Theorem 5.1 in [11]); but there are many other types of perturbations that could occur in CRNs in either an evolutionary setting, or in experimental/clinical settings. Given that many of the most complex signalling networks that arise in nature, such as signal transduction networks in mammalian cancer cells, will likely never be delineated in complete intricate detail in terms of elementary chemical reactions, the robustness-conferring mechanisms that self-assemble in biology require a completely new approach, an entirely new language, for understanding the fundamental structures of life's networks. In this paper we wish to offer just such an approach.
Let us look into this issue in more details. Suppose the dynamics is governed by a system of ODEs: where is the -dimensional state of the system and ( ) = ( & ( ), … , 0 ( )) are the rates of evolution of all the species. In order to exhibit an integrator, we need to find a real-valued function ( ) such the dynamics of ( ) = ( ( )) is given by the ODE: where ∇ is the gradient operator, • is the usual dot product and ( ) = ( )( 3 − ) is what the authors call the "RPA polynomial" for the RPA variable 3 . Suppose one is able to find nonlinear polynomials ℎ & , … , ℎ 0 (as in Theorem 1 of the paper) such that However this does not necessarily mean that a function ( ) can be constructed so that (0.1) holds. Observe that such a function ( ) would need to have the following partial derivatives Unless each ℎ 3 ( ) is a constant (in which case ( ) is linear), it is unclear why such a function ( ) would exist. Since the approach in the paper uses nonlinear elimination steps (via concatenating monomials), the existence of ( ) cannot be ascertained. In fact if for a pair , we have 6 ℎ % ( ) ≠ % ℎ 6 ( ) then certainly H(x) cannot be constructed with partial derivatives ℎ & , … , ℎ 0 . This discussion shows that while existence of a RPA polynomial may be necessary for integral control, it is not sufficient and hence this property does not characterise integral controllers (a main claim of the paper). This gap between necessity and sufficiency can be bridged if the authors can show that for RPA networks one always has 6 ℎ % ( ) = % ℎ 6 ( ) ∀ , .
Our response: Yes, the reviewer is entirely correct to point out that the nonlinear map we propose in this paper (Eq. 6) is not a coordinate transformation. But again, we don't actually seek a diffeomorphism ( ) of the type the reviewer describes, and we don't claim that our collection of elimination polynomials gives rise to a coordinate transformation (other than in the special case of constant elimination polynomials). We do apologise again for incorrectly referring to the nonlinear map as a nonlinear coordinate change in our previous rebuttal letter. This was sloppy (and incorrect) wording, and we have made absolutely certain that no such inaccuracies appear anywhere in our manuscript or in our SI.
It is true that the existence of an RPA polynomial, , contained in the ideal , 0 ]} is a necessary but insufficient condition for RPA. As our Theorem 1 makes clear, = ( 3 , 6 )( 3 − ) must generate the principal ideal in two variables, 7 ∩ ℝ[ 3 , 6 ], where 6 is (any) non-RPA-capable variable. This condition is both necessary and sufficient for RPA at 3 , under the assumption of stability. But nowhere do we claim that a CRN for which 8 3 , 6 : generates 7 ∩ ℝ[ 3 , 6 ] contains "an" integrator. This can only be true if the elimination polynomials are constants, as the reviewer points out. If non-constant elimination polynomials are required, our framework requires there to be multiple independent (linear) integrators, connected via concatenating monomials; the nonlinearity is thereby relegated entirely to the concatenating process, not to the identification of internal models. In general, the output for one internal model will be an input for another. This collection of linear coordinate maps is always obtainable from a decomposition of the special nonlinear transformation we identify. That is the viewpoint we develop here.
Again, we are re-interpreting the IMP from a completely different standpoint from the one considered in prior work, in which the existence of a nonlinear coordinate transformation, and an associated output-driven internal model, is normally at issue. Note that the RPA polynomial we reference in Eq. 6 (and Theorem 1 -the Two-Variable Kinetic Pairing Theorem) contains strictly two variables in contrast to (x) in (0.1) above, which involves an unspecified number of variables. We simply do not claim that the transformation that identifies the (two-variable) RPA polynomial is associated with a single integrator. Rather, we claim that this transformation is always able to obtain via a connected collection of linear integrators, each requiring a linear coordinate change, and each imposing RPA on some topologically important feature of the CRN.
Bin et al. in [1] offer a particularly insightful comment. They state in Section 3: 'The question that the IMP asks is, If a system … is seen experimentally to regulate against all inputs (in some class of functions, e.g. constants), then what can be said about its internal structure? Answers to this question may help guide experimentalists and modelers by ruling out putative mechanisms and suggesting a search for components responsible for adaptation'. The version of the IMP those authors go on to discuss in that paper (involving a coordinate transformation that identifies an output-driven internal model, within a feedback structure) provides one answer as to what can be said about the internal structure of RPA-capable networks. Our paper provides a different answer -one that reveals the universal properties, at the molecular level, of all RPA-capable networks.

Questions about the Kinetic Pairing result:
In the previous review round, I had proposed a possible RPA network and asked the authors how it satisfies the Kinetic Pairing result (Theorem 1). The network has two reactions given by and there may be several other reactions involving species ! , … , : that can be catalysed by 1 but they do not change 1. In the rebuttal letter, the authors consider certain instances of such a network, and claim that in these networks all possible steady states lie on the boundary on the positive orthant (e.g. ! or " is zero) and hence the networks are not RPA. In particular, the authors state in their rebuttal letter that: The two reactions proposed by the reviewer cannot impose RPA on X1 no matter which reactions are added to the original two. This is because the two given reactions are linearly independent, and the CRN formed from these reactions has deficiency zero. Moreover, the two complexes that contribute to the apparent RPA-polynomial in the reaction rate f1 are not members of the same terminal SCC. It is clear from our discussion of the decomposition of CRNs into algebraically independent subnetworks a decomposition which is closely tied to the partition of deficiency among independent subsets of reactions (see SI Section S1.3) -as well as our extensive discussions of rowspan polynomials (complex linear invariants) throughout the remainder of our supplement, that no possible addition of extra reactions can turn the original two into an RPA-conferring set. We briefly note in Section S4.4.1 (see also footnote in that section) the special case where a deficiency zero collection of reactions can engender RPA namely when the two relevant complexes reside in the same terminal SSC, and correspond topologically to a trivial (isolated) connector node. The pair of reactions suggested by the reviewer are clearly not of this type.
To test the authors' claim above, I tried to come up with instances of such networks (i.e. (0.2) + some reactions) which are RPA and yet do not have this issue of steadystate not being in the positive orthant. Please explain how the Kinetic Pairing result fits these networks.
Our response: Yes, in the previous round of review the reviewer provided us with the two reactions given in (0.2) above and, pointing out that such a reaction pair would produce a reaction rate , asked us how this could be reconciled to the claims of our Two-Variable Kinetic Pairing Theorem since this equation has the apparent form of an RPA polynomial, yet appears to have three rather than two variables. We pointed out that it is possible to have an RPA polynomial involving three (or more) species, but only two variables, in CRNs containing a 'boundary variable'. But more importantly, we went on to show that, for the particular reaction pair proposed by the reviewer equipped with the needed 'closure' reactions (for which we provided three different examples), the projection of the ideal generated by the collection of CRN rate equations ( 7 ) onto two variables (one RPA-capable, and one non-RPA-capable) produced the elimination ideal < ! " >. Hence, none of the three CRN examples could exhibit RPA.
Our Theorem 1 (the 'two-variable kinetic pairing theorem') states, simply, that a CRN is RPA-capable in some variable 3 exactly when the projection of the steady-state ideal, 7 , onto 3 and 6 (for any non-RPA-capable 6 ) is generated by a single polynomial of the form = 8 3 , 6 :( 3 − ). We call such polynomials 'RPA polynomials'. As we pointed out in our earlier rebuttal, it is not enough for a set of CRN equations to contain a reaction rate that has the form of an RPA polynomial (as is the case for all the reviewer's examples, both here and in the previous round of review, which all have We make some additional observations for each of the three new CRNs in turn below. As the reviewer notes, we had observed in our previous rebuttal that the two specific reactions originally proposed by the reviewer together have deficiency zero, and that the addition of the necessary closure reactions will not be able to increase this deficiency. As a consequence, that CRN could not exhibit RPA. In other words, the fact that those particular CRNs were not RPA capable was not tied to the specific choices of closure reactions that we provided as illustrative examples. Of course, if one begins with a different pair of reactions (i.e. not involving 9 and : ), as is the case for the first two new examples the reviewer provides, or adds in more than the needed closure reactions, as is the case for the third new example, then one can certainly increase deficiency, and thereby arrive at an RPA-capable set. But in this case, the RPA-capacity is tied inextricably to the altered structure arising from the different starting pair and/or the redundant closure reactions, with the attendant increase in deficiency associated with that altered structure, and is not due solely to the mere presence of the original pair (which appear to produce an apparent 'RPA polynomial' for 15 $ 12 in all cases). We explain this point in more detail in our ensuing analysis. We apologise if we unwittingly side-tracked the discussion by raising the issue of the original pair being characterized by a deficiency of zero, and that adding in the necessary 'closure' reactions to allow us to work with that original pair cannot increase deficiency. The point we were really trying to emphasize was that the mere presence of in the rowspan of the system, even for cases where ! " is a boundary variable, does not, in and of itself, imply RPA. Our Theorem 1 explains why this is the case.
Again, the fundamental issue in all six CRN examples (the three here, as well as the three from the previous round of review) is that the reaction rate for & does not actually generate the two-variable elimination ideal for these CRNs. In general, CRNs for which one of the reaction equations is an 'apparent' RPA polynomial may or may not exhibit RPA.
I start with a simpler version of this network with only three species (i.e. species 9 and : are absent) where I add inflow and outflow for species ! and " . So the overall network becomes This network has no boundary variables and so the variables are the same as species. It seems that this network will indeed show RPA for & with set-point & / ! . In Figure 1 I plot the simulated dynamics (rescaled to have the set-point of & as 1) for three randomly chosen values of initial states and rate constants, and one can see that RPA holds.
Our response: Here, the reviewer starts with a different version of the initial two reactions (rather than the two given in (0.2)), so this CRN can have a deficiency of one, even with only the minimum needed closure reactions. In any case, this CRN does indeed exhibit RPA because < >. We provide our Singular code and output below so that the reviewer can verify (see [1]): Whereas (0.2) includes two additional species, 9 and : , each assigned to a different reaction in the pair, thereby guaranteeing the linear independence of the reactions once the requisite closure reactions are added, this property is no longer present in the new reaction set. When organized correctly into its linkage classes, this CRN has six complexes, two linkage classes and a rank of three, resulting in a deficiency of one.
The CRN has just two non-terminal complexes, from which it follows that there is a polynomial of the form & ! " − ! & ! " in the rowspan of the system -a feature that was already obvious from the form of 15 $ 12 . But note that this fact, in and of itself, does not guarantee that RPA will obtain. The three examples considered in the previous round of review also contained 15 $ 12 = & ! " − ! & ! " but did not exhibit RPA. Note that the Shinar-Feinberg theorem presupposes that the system admits a positive steady-state; it does not, in and of itself, provide any method for checking that this will indeed be the case. Note also that for the cases previously considered in which ! " was a boundary variable, it was possible (in principle) for the generating RPA polynomial (referenced in our Theorem 1) to take the form The reason this particular CRN has a deficiency of one is that a linear combination of the reactions involving ! , " , 9 (in the second linkage class only) produces & → ∅, whereas a linear combination of the reactions involving & , " , 9 (involving the first linkage class, and thus involving the new pair of complexes ! + " and & ) produces ∅ → & . Thus, ! and " regulate both the production and the degradation of & . This makes the CRN a balancer module, where ! and " play the role of diverter species. Note that although the reaction involving : could be used in place of the reaction involving " in the argument above, these two reactions involve the same complexes (∅ and & ), so the 'redundancy' in the reactions does not result in an increase in deficiency (nor any distinct new cycles or feedforward actions). The same holds for using the reaction involving ; in place of the reaction involving 9 .
Next I modify this network to have the production of ! and " catalysed by & . So the new network becomes This network also appears to be RPA as shown by the simulations in Figure 2. Finally, I consider another variant of this network where we have species 9 and : that reversibly bind to each other to produce an inactive compound, and they catalyse the production of X2 and X3 respectively. Therefore the new network becomes This network also seems to be RPA, as shown in the simulations in Figure 3. → ∅, giving a deficiency of one, which creates the needed parallel routes for the production/degradation of & as discussed in the previous two examples. The reactions involving : , ; replicate the reactions involving " , 9 (respectively) using different sets of complexes, thereby increasing the deficiency from one to three. The reaction involving < replicates the reaction involving ' using the same set of complexes, and therefore does not contribute to any further deficiency increases. Our Singular code and output is provided below: This CRN has eleven complexes, four linkage classes, a rank of four, and therefore a deficiency of three. Once again, the added deficiency arises from the additional replication of the 'needed' closure reactions employing different sets of complexes.
These networks do not appear to satisfy the RPA characterisation given by this paper. Also consider Example 6.5.3 in [5] which does not seem to fit this result either. These counter-examples cast doubt on the correctness of the Kinetic Pairing Theorem on which the entire paper rests.
Our response: None of these networks constitute counter-examples. All of these networks satisfy the RPA characterization given by our paper. Our Theorem 1 (the 'Two-Variable Kinetic Pairing Theorem') correctly identifies that the three examples suggested here by the reviewer will exhibit RPA at & , since a two-variable RPA polynomial generates the principal elimination ideal 7 ∩ ℝ[ & , 6 ]. In addition, our Theorem 1 correctly identifies that the three similar examples considered in the previous round of review cannot exhibit RPA since the generator of the relevant elimination ideal is not an RPA polynomial (despite the fact that one of the rate equations has the apparent form of an RPA polynomial).
The reviewer mentions Example 6.5.3 in [5], which is the following: Although this example claims that 'the system shows ACR for both variables', this system does not describe an RPA-capable (or ACR-capable) CRN since there are no parameters, and no possible 'input' or disturbance, in this model. So there's really nothing for the network to adapt to here. Thus, there is no non-RPA-capable variable; both of the variables appear to exhibit 'RPA', since the system's steadystate is a single fixed point in ℝ ! (and again, that's not really what RPA is). The special type of RPA known as ACR normally pertains to CRNs where there are mass conservation relationships linking the molecular concentrations, such that the 'total' concentration of a particular type/class of molecule can be varied (by altering the initial conditions of the system, for instance). In any case, in order for a CRN to adapt to any kind of disturbance, there must always be at least one non-RPA-capable variable in the CRN (otherwise how can there possibly be any kind of 'internal model', or any network components that 'offset' the disturbance?) In any case, the definition of RPA adopted in our work (as general as this is!) does not accommodate 'non-adaptive' cases like this where there's simply nothing for the network to adapt to. We must respectfully point out that the author of [5] has been a little hasty in describing this as an ACR-capable network simply because there is a single positive solution & = 1, ! = 2.
In addition, the author of Example 6.5.3 claims that 'it is possible to find a reaction network modeled with mass-action kinetics, such that the associated system is The condition being referred to here ( 3 = 3 − 3 3 , with all coefficients of 3 , 3 non-negative) is the so-called Hungarian Lemma, which states (roughly) that for a polynomial dynamical system to be inducible by a CRN under the mass-action assumption, any term in a rate equation preceded by a negative sign must be divisible by the subject of the rate equation. In other words, for the & / equation, each term preceded by a negative sign must include & as a factor. But the Hungarian Lemma, strictly speaking, pertains to CRNs in which each chemical reaction occurs at a rate (normally noted as a kinetic parameter superposed on the associated reaction arrow), that is independent of the corresponding rates of all other reactions. Because Example 6.5.3 is not suitably parametrized in this sense, there is necessarily an ambiguity introduced in terms of which specific reactions might produce the various terms of the equations. In other words, there is not a one-to-one correspondence between reaction rates and putative CRN structure, as predicted by the Hungarian Lemma under the specified conditions. The rate & , for example, contains the term −2 & ! ; this could in principle From this point of view, we should not really even accept the two given equations as a valid CRN model. In any case, on both counts (not being ACR-capable, and not being a valid CRN model), this Example simply does not constitute any sort of counterexample to the claims of our paper.
We do very much appreciate the lengths to which the reviewer has gone in order to make absolutely certain that all our claims are 100% accurate and watertight. But we do hope the reviewer is now willing to accept the correctness of our assertions. Two variables is always the right number of variables for assessing RPA capacity, the mathematical reason for which is given in our proof to Theorem 1. For all RPAcapable CRNs, the geometric projection of their rate equations onto two variables is generated by a single polynomial of the form = 8 3 , 6 :( 3 − ). Any CRN that does not satisfy this property cannot exhibit RPA -no exceptions, un point c'est tout.
3. The set-point may not be a rational function of parameters: In many places in the main paper and the supplement, the system's set-point is said to be a rational function of biochemical parameters. While this holds for most examples, consider the following birth death network is not a rational function of c1 and c2? For more examples, see Examples 2.11 and 2.12 in [3]. Please revise as necessary.
Our response: We sincerely thank the reviewer once again for such a meticulous review of our work and for such careful attention to every possible technical detail. In this particular case, & ! is a boundary variable since & only appears in the CRN rate equations in the form of the monomial & ! . Therefore, the setpoint for this variable is & /2 ! , which is a rational function of parameters, as expected. With boundary variables accounted for in this way, we can be assured that the setpoint of any RPAcapable variable will be a rational function of parameters, since we are working with a commutative ring structure (multivariate polynomials) with coefficients taken over a field (real numbers, with all coefficients considered symbolically). We allow only addition and additive inverses (subtraction) and multiplication without inverses (i.e. no division) of the ring elements; the coefficients admit addition/subtraction as well as multiplication/division (due to the field structure). No radicals may be taken.
But we do completely appreciate the reviewer's point, and that the boundary variable technicality might be lost on some readers (despite the fact that we do highlight the concept quite prominently). The setpoint of the species concentration may indeed be an algebraic (rather than a rational) function of parameters in these special cases.
Examples 2.11 and 2.12 in [3] also contain a species setpoint that is an -th root of a rational function of parameters simply because the species in question contributes to an RPA-capable boundary variable in each case. In Example 2.11 -a continuation of Example 2.3 -( appears only in the form ( 0 . Likewise, in Example 2.12 -a continuation of Example 2.8 -$ appears only in the form $ 0 . We now add the following clarification of this matter to Remark 3 after the statement and proof to Theorem 1 (immediately following Definition 3, in Section S1.5 in our SI): We now also refer to this Remark in the main paper (at the end of the paragraph immediately following Eq. 6), where it is noted that is a rational function of parameters.

Questions on the Buchberger's Algorithm:
In the previous review round, I had raised the following issue regarding Buchberger's Algorithm The Gröbner basis algorithm to find the RPA polynomial may not terminate. It is mentioned that failure to terminate for a chemical reaction network (CRN) is a prima-facie evidence that the CRN does not exhibit RPA. However this is not mathematically shown. In any case, checking for nontermination of a method is impractical.
In response the authors state in their rebuttal letter that the termination is guaranteed because "all ideals of a polynomial ring (or any Noetherian ring, for that matter) are finitely generated -a result of central importance in algebraic geometry, formalized by the Hilbert Basis Theorem". While I agree that termination is guaranteed, I must point out that the reason for my confusion and for raising this termination issue is the following excerpt from the previous version of the Supplement (submitted in round 1): It seems that in this statement that authors meant "terminate in practical time". I understand that due to the "almost linear" nature of coordinate change, this algorithm would work well for RPA networks. However, does the final RPA polynomial produced by the method depend on the monomial ordering or other choices made by the method?
Our response: No, whether or not the elimination ideal referenced in our Theorem 1 is generated by an RPA polynomial is completely unrelated to the choice of monomial ordering. The form of the 'final RPA polynomial' is also unrelated to the choice of monomial ordering: it is the generator of a principal ideal. But it's important to bear in mind that only monomial orderings that can achieve the desired elimination are suitable for testing this property (i.e. an elimination ordering). One cannot apply some completely arbitrary monomial ordering and expect to obtain (any) elimination ideal. We do refer to the highly accessible text by Cox et al. [7] in both our main text and in our SI for the benefit of any reader who wants to understand these technicalities fully. More importantly, we provide sample code at the end of our SI, with a discussion of a particularly efficient elimination ordering (using command '(dp(n-2), dp (2))'see Section S5.1), should this be required for large and complicated CRNs.
Note that the reviewer's quotation from our original submission ('Failure of the Buchberger (or other Gröbner basis-computing) algorithm to terminate could thus be adduced as prima facie evidence that the CRN under consideration does not, in fact, have the capacity for RPA.') was removed from our revised submission, and is no longer included in our SI.
More importantly, if a network is not RPA, how many steps would the method need to confirm this non-RPA property? Such practical considerations must be discussed in the main text, and they are crucial for applying these ideas for characterising RPA in high-dimensional networks.
Our response: Of course there is no possible way to be specific about the number of steps that would be required to confirm the inability of some general (non-RPAcapable) CRN to exhibit RPA. If the algorithm is being applied completely mindlessly to a large and complicated CRN, with a random choice of two variables, then we quite agree that if the algorithm seems to be taking an inordinately long time to terminate, then the RPA-capacity or otherwise of the CRN is entirely unclear.
Now, the algorithm we present can be applied mindlessly, without any real understanding of the mathematical principles we discuss at length in our paper, and if the CRN under consideration is indeed RPA-capable, then such a completely mindless approach will be able to provide a confirmation of RPA capacity, along with an explicit identification of the setpoint. If the sole contribution of our paper were the development of an algorithmic test for RPA capacity, the reviewer would be quite correct to point out that this algorithm cannot easily distinguish between non-RPA-capable CRNs, and extremely large and complicated RPA-capable CRNs that could require a significant (and indeterminate) time-frame for the execution of the algorithm. But it's important to recognize that if the algorithm seems to be taking a long time to terminate, then our paper presents analytical approaches to check whether the CRN at hand has any potential for RPA capacity through a decomposition into independent subnetworks (wherever possible), and an analysis of deficiency in these subnetworks.
As an example, consider the CRN for the mammalian enzyme 6-phosphofructo-2kinase/fructose-2,6-bisphosphatase (PFK-2/FBPase-2), which operates bifunctionally to both activate and inactivate fructose-2,6-bisphosphate (F2,6BP), as discussed by Karp et al. [6] in terms of the potential RPA-capacity of F2,6BP. As we will show, this CRN cannot exhibit RPA. The CRN (presented in [6]), organized into its linkage classes, with & ≡ , ! ≡ -, " ≡ 6 , 9 ≡ --6 , : ≡ 2,6 , ; ≡ -2,6 , ' ≡ --2,6 , and < ≡ --F6P-2,6 , is given by Interestingly, Karp et al. [6] claim that the Gröbner basis implementation in Mathematica does not terminate for this CRN, although they give no information as to how they've ordered their variables, or which monomial ordering they chose (although they almost certainly would have used a lexicographic monomial ordering since, prior to our work, it was not known that an elimination ideal involving just two variables is all that's required to test RPA capacity; the lexicographic ordering is generally an extremely computationally expensive elimination ordering since it produces a full complement of elimination ideals, in vast excess of what's actually required to solve the RPA problem). We used Singular and used an efficient block monomial ordering (to project onto two variables) that allowed the algorithm to terminate in several seconds, from which we could confirm that this particular CRN is not RPA-capable (at : ≡ 2,6 or otherwise). We provide our Singular code below: As shown in our code, we have used an efficient block ordering here, with the two variables we wished to project onto in the first instance ( : and & ) in Block 1, and the remaining variables in Block 2, and with Block 2 ordered 'higher' than Block 1, and the efficient degree reverse lexicographic ordering ('dp') applied within each Block. (We describe the implementation of this efficient block ordering at the beginning of S5.1 in our Supplementary Information, where we indicate that the Singular command to implement this ordering is (dp(n-2), dp (2)), to be used in place of lp for lexicographic ordering). Only the beginning of the output is displayed above, as the first, second, and sixth polynomials in the Gröbner basis run to several pages in length due to the highly complicated coefficients for these polynomials. Upon careful scrutiny of the polynomials, it is evident that, for this CRN, 7  While this simple non-RPA-capable CRN in eight variables yields easily to the algorithm, it is certainly conceivable that vastly larger and more complicated CRNs could require a significantly longer timeframe for the algorithm to terminate. By way of general advice to the reader, if one allows the algorithm to run for several hours (e.g. overnight) using the most efficient implementation possible (as described above) and it still hasn't terminated after this period, then one should consider undertaking some additional analysis of the CRN (as we provide below) through (i) first decomposing into algebraically-independent subnetworks (if possible), and then (ii) considering where deficiency arises within these independent subnetworks, and whether the deficiency corresponds to the presence of parallel pathways and/or feedback cycles, as required for RPA. Consider also that one should have a reason to suspect that a CRN exhibits RPA. It is clear from the detailed principles we develop here, as well as our previous topological analysis of RPA-capable network architectures at the network macroscale [12], that RPA-capable networks are actually very 'special', with very specific structural requirements, and are therefore extremely rare in the space of all possible networks. If one has a completely arbitrary network, with no particular reason to suspect it might be able to exhibit this special type of robustness, then it almost certainly doesn't. In the specific case we analyse here, the CRN involves a bifunctional enzyme ( & ) which both upregulates and downregulates its target protein ( : ); a priori, this fact presents the possibility that this CRN might have the capacity for RPA (even if, in the final analysis, it actually doesn't). Now, recall that deficiency is a measure of the linear independence of the reactions of a CRN, relative to their distribution into the connected components (linkage classes) of a graph. In particular, the deficiency of a CRN is increased by one for every instance of a reaction being 'replicated' elsewhere in the network via a different set of complexes. Reversible reactions, for example, which duplicate a single reaction using the same pair of complexes, do not of themselves contribute to any deficiency increases. Recall also that, for an RPA-capable CRN, these replicated reactions must constitute either (i) a collection of parallel pathways, where the production/degradation (or interconversion between activation states) of the RPAmolecule must be orchestrated via different collections of reactions (involving different sets of complexes), or (ii) a cycle in the production/degradation (or interconversion between activation states) of a non-RPA-molecule (as regulated by the RPA-capable molecule), where the individual components of the cycle must be orchestrated via different collections of reactions (involving different sets of complexes). There's simply no other way for the right linear invariants to emerge! We highlight and discuss these principles repeatedly throughout the analysis of the worked examples in our SI.
Returning now to the non-RPA-capable CRN at hand: Note that the CRN for PFK-2/FBPase-2 has a deficiency of 5, comprising 14 complexes, 3 linkage classes and a rank of 6. By noting which reactions are 'replicated' (using different sets of complexes) to yield each of these 5 units of deficiency, it is straightforward to see that this CRN does not have the right structure to impose RPA on : (or any other molecule). First, we note that the eight species of the model are all intricately interconnected in the nineteen reactions of the CRN, such that no decomposition into independent subsets is possible (see discussion of this point in S1.3 of our Supplementary Information). Moreover, in assessing the linear independence of the CRN reactions, only one reaction in each pair of reversible reactions need be considered (for the reasons noted above pertaining to the deficiency-preserving nature of reversible reactions). Now, from a routine linear-algebraic analysis of the remaining reactions, it is easy to show that the reactions occurring at rates • Replace "consistutes" with "constitutes" on line 72 in the main text. Please run a spell check.
Our response: We thank the reviewer for spotting this typo, which has now been corrected. We have run a spell check, and there are no additional typos.
• The paper says that if a network is RPA the integrator is guaranteed to exist. For example the following on page 22 in the main text In principle, there should always exist some single nonlinear coordinate change to extract a single output-driven internal model (Fig. 9a) from systems rate equations, corresponding to a single integral of the systems tracking error (Fig. 9b) or the following on page 16 in the Supplement All classes of RPA, including ACR, thus require some form of integral control.
Please explain which version of IMP can be used to verify this existence. See [1] for a recent review on IMP.
Our response: Here we are referencing a long-standing idea in the systems biology literature (see, for example [8] and [9], and Section 3 of [1]) that, where RPA obtains, there 'should' exist some (generally nonlinear) coordinate change that recasts the system into integral feedback form (if necessary, if the system in question has a feedforward architecture, and is subjected to a disturbance at its 'diverter'). The integral in question should operate on the system's tracking error. Reference [8], in particular, provides a very detailed mathematical elaboration of these ideas, which hold for nonlinear systems under several technical assumptions (which are very mild in the context of CRNs). Both [1] and [8] apply these ideas to identify a nonlinear transformation that identifies an output-driven internal model for a simple CRN with a feedforward structure, through a recasting into feedback form. We reviewed this simple example in our response to the Reviewer's Point 1. Again, our paper provides a very different viewpoint, distinct from the notion that 'there should always exist some single nonlinear coordinate change to extract a single output-driven internal model … corresponding to a single integral'. Our reason for referencing the viewpoint of prior control theoretic literature [1,8,9] is to contrast it with our new completely different viewpoint, and not to present that prior viewpoint as axiomatic.
We now add in these three references [1,8,9] at the location noted by the reviewer, to emphasize the prior IMP viewpoint we draw upon here.
• On page 20 the authors state that "It is striking to note that the original form of the CRN (Fig. 8a) eludes the Shinar-Feinberg theorem, even though the CRN exhibits ACR and has a deficiency of one." However does it fit the results in [4]? Please comment on this. Also explain why the authors found the method for finding linear invariants "ad hoc" (lines 568-569 on page 25).
Our response: As noted in many of our previous responses, we suspect the reviewer is actually referring to Karp et al. [6], rather than Perez Millan et al. [4], in this comment. Certainly, once one knows which particular complexes (and hence, which monomials) one needs to include in any RPA/ACR-relevant invariants, one can certainly use the method developed Karp et al. [6] to compute those invariants. In this specific simple case, however, the invariants in question happen to be the actual rate equations / (the connector polynomial) and / (the balancer polynomial), so the method is hardly necessary.
Regarding our description of the Karp et al. [6] method for finding linear invariants as 'ad hoc', those authors themselves acknowledge that this method relies on knowing ahead of time which complexes (monomials) are relevant. As we had noted in one of our earlier responses, the authors specifically state that the 'examples (analysed in that paper) had already been analysed by other methods, so we had an idea of which invariants to expect and which subset of complexes to consider. For a new network such information may not be available, so how can non-trivial type 1 complex-linear invariants (simply, "invariants") be found?'. The authors go on to provide some practical suggestions for the reader on computing useful 'complex linear invariants', which may be helpful in studying new CRNs but certainly do not constitute a truly systematic procedure for rigorous application to general CRNs.
• On the Supplement page 15 it is stated that Mass-conservative CRNs therefore have no external stimuli or inputs, and can only be perturbed by altering the total abundances (or concentrations) of the constituent molecules -i.e. by altering the initial conditions.
Why cannot the perturbation come in the form of parameter variation, e.g. of a conversion reaction (that conserves mass).
Our response: Yes, the reviewer certainly does make an interesting point. In practice, conversion reactions are normally mediated by enzymes, so if the enzyme in question is not part of the network, we agree that the corresponding parameter (which reflects the concentration of the enzyme) could be perturbed. Almost any parameter of a CRN could be perturbed in principle. But the framework we develop here emphasizes that a CRN can adapt to perturbations to any parameter that does not appear in the RPA-variable's setpoint, regardless of whether the CRN is massconservative or not.